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SECTION  1 


INTRODUCTION  AND  SUMMARY 

1.1  Introduction 

This  report  has  been  prepared  as  the  Final  Report  of  a  preliminary 
study  of  the  possible  satellite  autonomous  navigation  accuracy  with  the 
SHAD  (Stellar  Horizon  Atmospheric  Dispersion)  concept  for  various  types  of 
representative  orbits. 

The  SHAD  concept  is  considered  to  be  one  of  the  most  promising  concepts 
yet  identified  for  achieving  high  accuracy  satellite  autonomous 
navigation.  Until  now,  the  orbits  of  most  satellites  have  been  determined 
by  using  ground  tracking  measurements  from  radar  systems,  telescopes,  etc. 
Unfortunately,  earthbound  tracking  suffers  from  many  limitations.  Loss  of 
a  tracking  station  during  wartime  could  prevent  ground  operators  from 
obtaining  the  necessary  orbit  determination  accuracy  for  certain  military 
satellites.  Ground  tracking  also  limits  civilian  satellites  since  the  high 
cost  of  tracking  strains  operating  budgets  and  contributes  to  the 
cancellation  of  many  otherwise  cost-effective  missions ^ ).  Furthermore, 
the  large  number  of  objects  in  orbit  taxes  the  capacity  of  existing 
monitoring  stations  and  delays  their  delivery  of  information  to  the 
operator  ^  ^ . 

In  the  near  future,  satellites  will  be  able  to  navigate  using  signals 
from  NAVSTAR  satellites^).  There  are  missions  being  planned,  however, 
where  even  this  degree  of  self-sufficiency  may  be  inadequate ^ 4 ) . 


To  meet  these  needs,  satellites  can  be  built  which  estimate  their 


positions  without  recourse  to  either  ground  stations  or  other  satellites. 


Such  an  "autonomous  navigation  system"  would  allow  a  satellite  to  carry  out 


many  functions  without  instructions  from  the  earth.  Many  people  have 


suggested  ways  to  navigate  satellites  autonomously;  a  few  of  these  methods 


have  already  been  implemented.  However,  the  proposed  and  existing 


navigators  have  not  been  widely  accepted,  in  part  because  of  their  inherent 


complexity,  poor  accuracy,  or  exorbitant  computational  requirements  1 b > . 


The  SHAD  concept  takes  advantage  of  the  atmosphere's  optical 


properties.  Viewed  from  orbit,  a  star  passing  behind  the  earth's  upper 


atmosphere  will  appear  to  be  shifted  upward  from  its  true  position.  This 


refraction  is  a  property  of  starlight  frequency  and  atmospheric  density. 


If  a  refraction  measurement  were  to  be  made  on  a  known  star  near  the 


horizon,  this  would  provide  information  on  the  direction  of  that  portion  of 


the  earth's  horizon  with  respect  to  the  satellite  in  inertial  coordinates. 


There  are  essentially  two  ways  of  measuring  starlight  refraction.  One  is 


is  to  measure  the  refraction  directly  such  as  with  a  star  tracker.  The 


other  is  to  measure  the  difference  in  refraction  (i.e.  the  dispersion) 


between  two  wavelengths  of  the  starlight  such  as  red  and  blue.  The  latter 


approach  is  utilized  by  the  SHAD  sensor  and  there  are  several  advantages  to 


performing  the  measurement  in  this  manner.  However,  it  should  be  noted 


that  the  above  two  methods  are  merely  different  ways  of  measuring  the  same 


basic  phenomenon.  The  relationship  between  refraction  and  dispersion  is 


very  well  known  for  air.  Consequently,  accurate  measurements  of  one 


automatically  imply  accurate  knowledge  of  the  other. 


t1.  V. V. 
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Use  of  refraction  or  dispersion  measurements  for  satellite  navigation 


will  require  access  to  an  accurate  atmospheric  refraction  model  and  a  star 
catalog*  Each  measurement  will  provide  information  about  a  particular 
component  of  satellite  position.  By  performing  these  measurements  on 
different  stars  throughout  the  orbit,  a  complete  determination  of  the 
satellite  position  and  velocity  can  be  made.  The  accuracy  of  this 
determination  will  depend  upon  many  paramenters  such  as  measurement  error, 
number  and  direction  of  star  sightings,  and  type  of  satellite  orbit.  One 
of  the  primary  objectives  of  this  study  is  to  provide  some  indication  of 
the  possible  navigation  accuracy  and  to  indicate  the  sensitivity  of 
performance  to  variation  in  the  most  important  parameters. 


1 .2  Background 

During  the  Apollo  program,  the  Charles  Stark  Draper  Laboratory  (then 
the  M.I.T.  Instrumentation  Laboratory)  investigated  many  techniques  for 


performing  orbit  navigation  in  the  Command  and  Lunar  Modules’ 


Early 


in  this  effort,  mission  planners  felt  that  ground  tracking  was  not 
sufficiently  accurate  or  reliable  enough  to  permit  navigation  of  manned 
spacecraft  over  large  distances.  What  was  desired  was  a  system  that 
allowed  the  crew  to  determine  its  position  without  depending  upon 
earth-based  systems.  Several  navigation  schemes  considered  for  this 


purpose  involved  the  sensing  of  the  Earth's  horizon  either  by  measuring  the 


horizon  radiance  in  some  portion  of  the  electromagnetic  spectrum  (e.g.  IR, 
UV,  visible,  etc.)  or  by  measuring  the  effect  on  starlight  (i.e. 
attenuation,  refraction,  etc.)  when' it  grazes  the  atmosphere ^ 7 ) . 


Although  some  of  these  approaches  were  very  attractive,  the  technique 
finally  adopted  was  to  use  a  sextant-like  device  which  measures  the  angle 
between  an  earth  or  lunar  feature  (e.g.  landmark,  limb,  etc.)  and  a  known 
star.  By  the  time  Apollo  was  launched,  however,  ground  tracking  had 
improved  to  the  point  where  the  autonomous  navigation  system  could  be 
relegated  to  back-up  status. 

In  1975,  interest  in  horizon  sensing  by  starlight  refraction  was 
renewed  at  CSDL  in  a  Navy  contract  concerned  with  investigating  all  known 
techniques  for  autonomous  navigation  of  vehicles  traveling  through 
space Though  many  people  expressed  doubts  about  accurately  modeling 
and  measuring  starlight  refraction,  a  few  researchers  felt  that  the  method 
had  potential  since  better  star  trackers  were  now  available  for  this 
application  and  there  was  no  indication  that  the  refraction  modeling  was 
serious.  Several  sensor  designs  were  conceived  for  accurately  measuring 
starlight  refraction;  and  a  somewhat  different  type  of  sensor  was  jointly 
proposed  by  ONR,  CSDL  and  Un.  of  MD  in  1979  to  measure  the  starlight 
dispersion  instead  of  the  refraction.  Although  dispersion  is  simply  the 
variation  in  refraction  with  wavelength,  previous  development  and  tests  of 
a  two-color  ref ractometer  (TCR)  at  the  Un.  of  MD  indicated  that  this  device 
could  measure  the  dispersion  between  red  and  blue  light  with  very  high 
accuracy.  For  autonomous  navigation,  this  device,  with  appropriate 
modifications  for  star  acquisition,  tracking,  etc.,  is  denoted  as  the  SHAD 
sensor. 

i 

Having  Jsdent^g^d  several  ways  of  accurately  measuring  refraction  or 
dispersion,  there  was  still  the  problem  of  improving  our  knowledge  of 
refraction  modeling  and  the  general  behavior  of  the  basic  phenomenon.  A 


significant  step  was  achieved  in  this  direction  when  CSDL  conducted  a 
survey  of  existing  satellites  in  late  1979  to  determine  if  any  could 
provide  real  data  on  stellar  atmospheric  refraction.  Two  NASA  satellites 
(0A0-3  and  HEAO-2)  were  found  to  have  this  capability  and  CSDL  was  granted 
guest  observer  status  to  obtain  this  data.  Successful  observations  of 
refraction  were  obtained  from  OAO-3  (Orbiting  Astronomical  Observatory)  on 
May  14  and  July  11,  1980.  Further  experimentation  with  OAO-3  was  not 
conducted  since  it  was  found  that  HEAO-2  (High  Energy  Astronomy 
Observatory)  had  inadvertently  been  making  such  observations  with  high 
accuracy  since  it  was  launched  in  November  1978.  It  is  estimated  that 
several  thousand  observations  were  made  by  HEAO-2  before  it  was  deactivated 
in  April  1981.  Atleast  700  of  these  observations  were  retrieved  and 
analyzed  to  some  extent  by  CSDL,  and  a  limited  number  of  cases  were 
processed  with  the  highest  available  accuracy,  using  more  accurate  NASA 
determinations  of  the  satellite  orbit  with  ground  tracking  data.  The 
results  of  that  effort  indicate  that  the  concept  has  a  very  high  navigation 
accuracy  potential.  The  stratospheric  density  variations  indicated  by  the 
more  accurately  processed  data  were  much  less  than  those  indicated  by 
meteorological  balloon  and  rocket  observations. 

1 . 3  Summary 

The  primary  objective  of  this  study  was  to  generate  some  indication  of 

the  autonomous  navigation  accuracy  that  can  be  achieved  with  the  SHAD 

« 

concept  for  various  representative  satellite  orbits.  This  was  accomplished 
by  conducting  a  covariance  analysis  utilizing  a  Kalman  filter  to 


sequentially  process  the  measurements  as  they  occur  during  orbit 

operation.  ^  Past  experience  has  indicated  that  a  significant  amount  of 

•  \ 

qualitative  and  reasonably  accurate  performance  data  can  be  conveniently 
generated  with  this  type  of  simulation.  Such  data  not  only  provides  a 
fairly  good  indication  of  performance  but  also  minimizes  the  number  of 
cases  that  one  may  later  wish  to  analyze  more  accurately  with  a  full-state 
simulation  utilizing  elaborate  gravity  and  atmospheric  models,  etc. 

The  performance  results  obtained  in  this  investigation  indicate  that 
high-accuracy  autonomous  navigation  is  possible  using  relatively  few 
stellar  refraction/dispersion  measurements  per  orbit.  In  the  baseline 
cases  examined,  steady  state  position  error  standard  deviations  were  often 
much  less  than  100  meters  which  is  superior  to  most  satellite  navigation 
systems.  Initial  large  uncertainties  in  the  navigator's  knowledge  of 
satellite  position  and  velocity  were  rapidly  overcome  within  the  first 
orbit  of  navigation  and,  thereafter,  the  performance  slowly  converged  to 
steady-state  levels  in  about  five  or  six  orbits. 

Increasing  the  number  of  stars  sighted  per  orbit  improves  overall 
accuracy  but  near-optimal  performance  can  be  achieved  with  relatively  few 
sightings.  Widening  the  azimuth  span  is  valuable  for  optimizing  accuracy 
in  the  normal  direction.  However,  for  a  given  number  of  visible  stars  in 
the  sky,  the  best  performance  is  obtained  by  making  star  sightings  both 
fore  and  aft  of  the  satellite's  direction  of  motion. 

The  high  accuracies  obtained  are,  of  course,  due  in  large  part  to  the 

i 

low  measurement  variance  used  in  the  simulations.  Preliminary  refraction 
data  from  the  satellite  HEAO-2  show  that  the  atmosphere  is  sufficiently 
stable  in  the  tropics  to  warrant  this  low  uncertainty.  As  more  is  learned 
about  the  upper  atmosphere,  this  figure  will  certainly  be  revised  to  better 


studies. 

There  is  some  indication,  based  on  other  independent  navigation 
studies,  that  the  process  noise  levels  used  in  the  present  study  are  less 
than  what  they  should  have  been  for  the  relatively  simple  gravity  modeling 
assumptions  made  for  the  navigator.  If  this  is  the  case  then  the 
performance  results  in  this,  report  are  somewhat  optimistic  for  a  simple 
navigator.  However,  they  would  be  more  correct  for  a  navigator  utilizing 
more  complete  models  of  gravity,  drag,  etc  which  is  not  only  possible,  but 
probably  desirable. 

As  previously  noted  the  present  study  was  essentially  a  covariance 
analysis  which  provides  a  reasonably  good  indication  of  performance  but 
does  not  provide  the  accuracy  of  a  full-state  simulation.  When  one 
considers  the  rather  complex  environment  this  concept  must  operate  in,  it 
would  seem  most  appropriate  to  develop  a  full-state  simulation  in  any 
future  study  of  this  concept. 
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SECTION  2 


SATELLITE  NAVIGATION  BY  STELLAR  REFRACTION 

2.1  Atmospheric  Refraction 

The  passage  of  starlight  through  the  earth's  atmosphere  bends  the  rays 
inward  (see  Figure  2-1).  Viewed  from  orbit,  the  star's  apparent  position 
remains  on  the  horizon  long  after  it's  true  position  has  "set".  In  Figure 
2-1  the  refracted  ray  observed  by  the  satellite  appears  to  be  grazing  the 
earth  at  an  apparent  height  ha  but  actually  grazes  the  earth  at  a 
slightly  lower  height  denoted  as  hg. 

Refraction  is  strongest  near  the  surface,  becoming  progressively  weake 
at  higher  altitudes.  In  effect,  the  atmosphere  acts  like  a  prism,  refract 
ing  and  dispersing  the  starlight  passing  through  it.  A  ray  of  starlight 
passing  through  the  spherical  shell  of  the  atmosphere  encounters  the 
gradient  in  air  density  which  determines  the  amount  by  which  the  starlight 
is  bent.  Densities  near  the  Earth's  surface  are  known  to  be  closely 
described  by  an  exponential  function  of  altitude  such  as  the  following: 


/  h'ho\ 

=  Po  exp  \-  — - ) 


(2-1  ) 


where  p  is  the  density  at  altitude  h,  pQ  is  the  density  at  some  other 
altitude  hQ,  and  H  is  the  density  scale  height  which  is  defined  by 


H  = 


R  T 

m 


M  g  +  R  /dT  \ 

°  g  / _ E  1 

\dh  / 


(2-2) 


where  Rg  is  the  universal  gas  constant,  Tm  is  the  molecular  scale 
temperature,  M0  is  the  sea-level  value  of  the  molecular  weight  of  air. 


and  g  is  the  acceleration  due  to  gravity.  Another  expression  for  the 


density  scale  height  which  may  be  easier  to  visualize  is  obtained  by 


inverting  Equation  2-1  to  obtain: 


H  =  (h-hc)/  m(p0/p) 


(2-3) 


Since  the  gravitational  acceleration  g  and  the  molecular  scale  tempera 


ture  Tm  change  with  altitude  for  a  realistic  atmosphere,  one  can  expect 


some  change  in  H  with  altitude.  However,  the  change  in  H  over  a  limited 


altitude  range  is  generally  not  very  large. 


If  the  value  of  H  at  some  altitude  hg  is  assumed  to  be  the  same  at 


all  higher  altitudes,  an  approximate  value  of  the  starlight  refraction 


angle  (7»8)  is  given  as  follows: 


R  -  (yg-1) 


2x  ,r  +  h 


(2-4) 


where 


=  total  refraction  (rad) 


=  tangent  (grazing)  height  of  refracted  ray  (km) 


=  earth  radius  (km) 


index  of  refraction  at  height  hg 


=  density  scale  height  at  hg  (km) 


According  to  Gladstone  and  Dale's  law,  the  index  of  refraction  (y  )  is 


related  to  the  air  density  (p)  as  follows: 


y  -  1  =  k  (X  )p 


(2-5) 


where  k(X  )  is  denoted  as  the  dispersion  parameter  and  is  a  function  of  only 


the  wavelength  of  the  light.  Therefore,  the  refraction  angle  in  Equation 
2-4  may  be  expressed  in  terms  of  the  air  density  (pg)  at  the  grazing 
height  (hg)  instead  of  the  index  of  refraction  (yg): 


R 


k(X)p 


2it  (r  + 


h  ) 

_2_ 


H 


(2-6) 


where  an  expression  for  k(X)  may  be  obtained  by  considering  Equation  2-5 
for  air  at  standard  temperature  and  pressure.  According  to  Edlen^9^  the 
index  of  refraction  (ps)  for  air  at  standard  temperature  and  pressure  is 
related  to  light  wavelength  as  follows: 


<wa-i  >  io8 


6432.8  ♦  1949810  +  25540 

146  -  1/x2  41  -  l/^2 


(2-7) 


where  X  is  the  wavelength  in  microns.  The  air  density  at  standard  tempera¬ 
ture  and  pressure  is  1225  g  m-^.  Substituting  the  above  values  for  air 
density  and  index  of  refraction  into  Equation  2-5  and  solving  for  k ( X ) 
yields 


k(X  ) 


1  x  10 


-8 


5.2513  + 


2408 


20.849 


(2-8) 


146  -  1/ 


41  -  1/ 


To  illustrate  how  well  Equation  2-6  approximates  the  refraction  of  a 

t 

realistic  atmosphere  such  as  the  1976  U.S.  Standard  Atmosphere 1 0 ^  the 
results  with  Equation  2-6  are  compared  in  Table  2-1  with  those  generated  by 
accurate  ray  tracing^11 K  The  refraction  angle  given  by  Equation  2-6  for 
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a 


a 
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Table  2-1.  Refraction  angles  for  the  1976  U.S.  Standard  Atmosphere 


Refraction  Angle  (arcsec)  Computed 

At  X  =0.7  microns  By: 

Grazing 

Height 

(km) 

Density 

Scale 

Height  (km) 

Equation  2-6 

Accurate  Ray  Trace 

20 

6.225 

331  .7 

333.8 

25 

6.366 

147.9 

148.1 

30 

6.519 

67.2 

•67.3 

35 

6.508 

30.9 

30.7 

40 

6.897 

14.2 

14.1 

45 

7.288 

6.8 

6.7 

50 

8.047 

3.4 

3.3 

-  19  - 
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each  grazing  height  was  computed  using  the  density  and  density  scale  height 
indicated  for  that  grazing  height  in  the  U.S.  Standard  Atmosphere.  Note  in 
Table  2-1  that  the  density  scale  height  of  the  standard  atmosphere  does 
vary  slightly. 


2.2  Atmospheric  Dispersion 

Atmospheric  dispersion  represents  the  natural  variation  in  atmospheric 
refraction  with  wavelength  in  accordance  with  the  dispersion  parameter  k(x  ) 
given  in  Equation  2-8.  The  passage  of  starlight  throuc^i  the  earth's  atmo¬ 
sphere  will  cause  the  blue  rays  to  bend  more  than  the  red  rays  (see  Figure 
2-2).  A  measure  of  the  dispersion  is  provided  by  measuring  the  difference 
in  refraction  between  the  red  and  blue  rays.  Since  the  optical  properties 
of  natural  air  are  well  known,  there  is  a  well  known  relationship  between 
dispersion  and  refraction  for  two  wavelengths,  and  one  may  convert  either 
type  of  measurement  to  its  counterpart  without  significant  loss  of 
accuracy. 

For  a  satellite  at  1000  km  altitude  the  values  of  dispersion  observed 
between  wavelengths  of  0.35  and  0.7  microns  for  different  grazing  heights 
in  the  1976  U.S.  Standard  Atmosphere  are  shown  in  Table  Note 

that  the  dispersion  values  are  small  in  comparison  to  the  refraction 
values.  However,  it  should  be  noted  that  a  two-color  ref ractometer  has 
been  developed  which  can  very  accurately  measure  the  dispersion,  and  also 
has  other  design  features  of  merit  in  the  intended  application.  A  detailed 
treatment  of  dispersion  and  its  relationship  to  refraction  is  given  in 
Section  2.4.6. 


**  a'  T 


Table  2-2.  Approximate  Values  of  Stellar  Refraction  and  Dispersion  at 
Wavelengths  of  0.35  and  0.7  Microns  for  the  1976  U.S. 
Standard  Atmosphere  and  a  Satellite  Altitude  of  1000 
Ki lometers 


Grazing 

Height 

(km) 

Refraction 
at  0.7  microns 
(arcsec ) 

Dispersion  Between 

0.35  and  0.7  microns 
(arcsec ) 

20 

333.76 

6.44 

25 

148.13 

3.88 

30 

67.27 

2.12 

35 

30.69 

1  .05 

40 

14.07 

0.50 

45 

6.70 

0.25 

50 

3.34 

0.12 
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2.3  Improved  Atmospheric  Density  Models 

The  accuracy  of  horizon  sensing  with  the  SHAD  concept  will  depend  on 
our  knowledge  of  the  stratospheric  density  profile.  Until  now,  most  of  the 
data  on  stratospheric  density  has  been  obtained  indirectly  through  measure¬ 
ments  of  temperature  and  pressure  made  by  instruments  carried  aloft  in 
high-altitude  balloons  and  rockets.  Extensive  temperature  data  has  also 
been  obtained  with  vertical  sounders  on  satellites.  The  accuracy  of  the 
above  measurements  has  generally  been  around  2  to  3  percent  and,  conse¬ 
quently,  places  a  limit  on  how  well  the  density  can  be  estimated  and 
modeled. 

The  variability  of  stratospheric  density  is  known  to  be  a  function  of 
at  least  latitude  and  season.  Stratospheric  observations  with  meteorologi¬ 
cal  rockets  and  balloons  indicate  a  density  uncertainty  at  25  km  altitude 
of  about  1.3  percent  in  the  tropical  region ^ 1 2» 1 3  ) .  For  the  summer  hemi¬ 
sphere  the  indicated  density  uncertainty  may  increase  by  a  factor  of  two  or 
three  when  going  from  the  equator  to  the  pole.  However,  the  winter  hemi¬ 
sphere,  especially  in  the  north,  may  experience  local  variations  as  high  as 
10  percent.  It  should  be  noted  that  most  of  the  above  data  is  based  on 
temperature  measurements  which  have  been  combined  with  pressure  estimates 
to  obtain  density  estimates.  The  errors  in  the  density  estimates  for  a 
given  observation  have  generally  been  no  better  than  one  or  two  percent. 

More  accurate  stratospheric  density  measurements  on  a  worldwide  basis 
would  disclose  certain  systematic  changes  in  the  stratosphere,  heretofore 
not  known  or  firmly  established  because  of  limitations  in  existing  data. 
Such  an  endeavor  is  being  proposed  as  one  of  the  objectives  of  a  SHAD 


satellite  experiment  where  either  dispersion  or  refraction  measurements 


would  be  used  in  combination  with  very  accurate  independent  knowledge  of 


satellite  position  to  estimate  the  observed  density.  Inspection  of 


Equation  2-6  indicates  that  refraction  is  essentially  a  direct  indication 


of  density  at  the  grazing  altitude.  Although  the  derivation  of  density 


from  refraction  measurements  is  not  quite  that  straight-forward,  relatively 


simple  techniques  have  been  developed  and  are  new  being  used  to  estimate 


density  from  actual  observations  of  refraction  made  by  the  satellite 


HEAO-2.  Although  this  satellite  was  never  designed  to  make  stellar 


atmospheric  refraction  observations,  CSDL  discovered  in  late  1979  that 


HEAO-2  had  inadvertently  been  making  such  observations  with  its  two 


guide-star  trackers  since  its  launch  date  (November  1978),  sometimes  as 


often  as  once  per  orbit. 


As  of  this  writing,  CSDL  has  processed  about  700  HEAO-2  observations 


taken  throughout  the  world  between  +  45  degrees  latitude.  Detailed  analy¬ 


sis  of  about  140  observations  taken  in  the  tropical  region  at  different 


times  of  the  year  indicate  that  the  stratospheric  density  variability  is 


less  than  1  percent.  The  accuracy  of  the  HEAO-2  data  in  providing 


estimates  of  density  is  primarily  limited  by  the  accuracy  of  the  HEAO-2 


orbit  determination  from  ground  tracking  data. 


2.4  Naviqation  by  Stellar  Refraction  or  Dispersion 


satellite  navigation.  This  is  then  followed  by  a  discussion  of  dispersion 


and  how  it  may  be  converted  to  an  equivalent  indication  of  refraction  for 
the  purposes  of  satellite  navigation. 


2.4.1  Basic  Navigation  Concept 

Measurements  of  atmospheric  refraction  or  dispersion  for  stars  near  the 

earth's  horizon  provide  an  indication  of  the  direction  of  the  horizon  with 

respect  to  inertial  space  which,  in  turn,  provides  an  indication  of  the 

satellite  position  with  respect  to  an  earth-centered  inertial  coordinate 

system.  This  is  graphically  illustrated  in  Figure  2-3  for  the  case  where  a 

satellite  observes  a  particular  value  of  refraction  (e.g.  150  arcsec)  for  a 

known  star.  If  one  assumes  a  spherically  symmetric  atmosphere,  all  of  the 

starlight  refracted  by  a  given  amount  will  define  a  conical  surface 

extending  out  into  space  and  having  an  axis  passing  through  the  center  of 

the  earth  in  the  direction  of  the  star.  Observation  of  this  particular 

value  of  refraction  by  a  satellite  indicates  that  it  is  somewhere  on  the 

surface  of  the  cone.  By  repeating  the  same  type  of  observations  on  stars 

in  different  directions,  the  navigator  can  determine  the  satellite's 

complete  position  by  essentially  solving  for  the  intersections  of  the 

various  cones.  It  should  be  noted,  however,  that  the  complexity  of  solving 

for  cone  intersections  is  seldom  necessary  since  the  navigator  usually  has 

sufficiently  accurate  knowledge  of  the  satellite  position  before  each 

measurement  to  permit  it  to  use  a  much  simpler  technique  for  position 

« 

update.  At  the  time  of  measurement,  the  navigator  will  usually  have  a 
prior  estimate  of  the  satellite  position  which  will  be  in  the  vicinity  of 
some  small  region  of  the  cone.  Since  the  measurement  indicates  that  the 


STARLIGHT  APPARENT  HEIGHT 


satellite  is  on  the  cone,  the  most  probable  position  is  that  point  on  the 


cone  closest  to  the  estimated  position.  Thus,  the  navigator  should  update 
along  the  perpendicular  from  the  estimated  satellite  position  to  the  cone's 
surface. 

It  is  apparent  that  the  above  update  provides  positional  information  in 
only  one  dimension.  However,  similar  updates  for  horizon  stars  in  other 
directions  throughout  the  orbit  will  provide  a  complete  update  of  position 
and  velocity. 

2.4.2  Refraction  Cone  Dimensions 

The  actual  dimensions  of  the  refraction  cone  are  a  function  of  the 
particular  value  of  refraction  being  considered  and  the  atmospheric  density 
model  assumed.  For  every  star  there  is  an  infinite  set  of  refraction  cones 
corresponding  to  an  infinite  set  of  refraction  angles  and  corresponding 
grazing  altitudes.  The  apex  angle  of  a  given  cone  is  equal  to  twice  the 
refraction  angle  associated  with  that  cone,  while  the  base  dimension  is  a 
function  of  the  earth's  diameter  and  the  grazing  height  of  the  refracted 
light.  The  geometry  associated  with  the  base  of  the  cone  is  shown  in 
Figure  2-4.  Here  it  is  seen  that  the  base  radius  might  best  be  represented 
as  a  function  of  the  apparent  ray  height  (ha )  instead  of  the  grazing  ray 
height  (hg).  For  triangle  AOB,  the  base  radius  (b)  is: 

b  =  (re  +  ha)/cos(R)  (2-9) 

* 

Also  note  in  Figure  2-4  that  the  base  radius  may  be  expressed  as 

b  =  r  +  h  +  a  (2-10) 

e  a 
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Figure  2-4.  Cone  base  geometry 


where  the  term  (a)  is  a  small  distance  which  can  be  shown  to  be  the 
following  by  equating  Equation  2-9  to  Equation  2-10: 

a  =  (  - 7TT-  -  l)  (r  +  h  ) 

'  cos(R)  1  ^  e  a' 

For  a  typical  refraction  angle  of  150  arc-seconds  at  ha  =  25  km,  the  term 
(a)  will  have  a  value  of  1.69  meters  which  can  usually  be  neglected  so  that 

b  -  r  +  h  ( 2-1 1 ) 

e  a 

It  is  important  to  note  that  the  distance  term  (a)  and  the  refraction  (R) 
have  been  greatly  enlarged  in  Figure  2-4  for  ease  of  illustration.  The 
refraction  values  of  interest  in  this  application  are  actually  quite  small 
and  one  could  essentially  visualize  the  refraction  cone  as  being  a 
cylinder . 

Having  expressed  the  radius  of  the  cone's  base  as  a  simple  function  of 
ha  (Equation  2-11),  there  is  the  question  of  how  ha  relates  to  hg  and 
the  associated  refraction  for  a  given  atmosphere. 

If  the  atmosphere  is  spherically  stratified  in  the  region  where  refrac¬ 
tion  occurs,  application  of  Snell's  law  leads  to  the  following  important 
relationship  for  any  point  on  the  ray  path: 

y  r  sin(Z)  =  Constant  ’  (2-12) 

where  y  is  the  index  of  refraction  at  a  given  point,  r  is  the  radial 
distance  of  that  point  from  the  center  of  the  earth,  and  Z  is  the  zenith 


angle  of  the  ray  at  that  point.  The  constant  is  sometimes  referred  to  as 
the  impact  parameter  and  may  be  determined  by  applying  Equation  2-12  to  the 
grazing  point  G  of  Figure  2-4  where  Z  =  90°,  r  =  re  +  hg,  and  y  =  yg, 

**e*'  Constant  =  y  fr  +  h  )  (2-13) 

Mg  i  e  gJ 


Substituting  the  above  value  of  the  constant  into  Equation  2-12  and  apply¬ 
ing  the  equation  to  the  satellite  point  S  of  Figure  2-4,  where  Z  =  Zs  and 
y =1 ,  yields  the  following: 

sin(Z  )=y  fr  +  h  )  /  r  (2-14) 

s  g  '  e  g'  s 

Also  note  for  triangle  SAO  of  Figure  2-4  that 

sin(Z  )  =  fr  +  h  ]  /r  (2-15) 

s  e  a}  s 


Therefore 


(2-16) 


or 


g  g 


(2-17) 


Note  that  the  above  relationship  between  ha  and  hg  requires  only  the 
index  of  refraction  at  hg.  Also,  note  that  this  relatively  simple  rela¬ 
tionship  applies  to  any  spherically  stratified  atmosphere,  whether  it 
varies  exponentially  with  altitude  or  otherwise.  Since  yg  can  be 


expressed  in  terms  of  pg  and  k(A)  in  accordance  with  Gladstone  and  Dale's 
law.  Equation  2-17  may  also  be  expressed  as  follows: 


h  =  k{\)  p  (r  +  h  )  +  h 
a  Kg'  e  g; 


(2-18) 


At  the  altitudes  of  interest  in  the  present  study  (i.e.  20  km  and  greater), 


the  first  term  on  the  right  of  Equation  2-18  will  be  very  small,  causing 


ha  to  be  slightly  larger  than  hg.  An  indication  of  the  relative  magni¬ 
tudes  of  ha  and  hg  is  given  in  Table  2-3  for  two  different  representa¬ 


tive  atmospheres  (1976  U.S.  Standard  and  1962  Tropical).  Note  that  the 


values  of  ha  obtained  with  Equation  2-18  are  the  same  as  those  obtained 


by  accurate  ray  tracing.  This  is  to  be  expected  since  the  ray-trace 


results  were  for  two  atmospheres  which  were  assumed  to  be  spherically 


stratified.  Also  note  that  the  relationship  between  ha  and  ha  is 


essentially  the  same  for  both  of  the  atmospheres  considered  in  Table  2-3, 


even  though  there  are  larger  proportional  differences  in  the  corresponding 


densities  and  refraction  angles.  At  a  grazing  height  of  20  km  it  is  seen 


that  the  density  is  about  7  percent  higher  for  the  tropical  atmosphere. 


However,  the  corresponding  apparent  height  is  only  9  meters  higher. 


2.4.3  Relationship  Between  Satellite  Position  and  Refraction  Cone 


Having  established  the  dimensions  of  the  refraction  cone  as  a  function 


of  refraction  angle  and  apparent  ray  height,  the  question  is  now  one  of  how 


a  refraction  measurement  can  be  used  to  relate  the  satellite  position  to 


the  cone.  This  problem  is  facilitated  by  introducing  a  new  quantity  re¬ 


ferred  to  as  the  vacuum  tangent  height  (hv)  of  the  star  LOS  (i.e.  the 


tangent  height  of  the  star  LOS  if  no  atmosphere  were  present).  It  will  be 


shown  that  the  vacuum  tangent  height  provides  the  "connecting  link"  between 


the  satellite  position  and  the  refraction  cone  since  it  may  be  expressed 


-  31  - 


Table  2-3.  Refraction  data  at  X  =  0.7  microns  for  the  1976  U.S 


Standard  and  1962  Tropical  Atmospheres 


hg 

(km) 

ha  (km)  Estimated  By 

Density  (pg) 

(g/m3  ) 

R 

(arcsec ) 

Equation  2-18 

Ray  Tracing 

U.S.  Standard: 

20 

20.128 

20.128 

88.910 

333.76 

25 

25.058 

25.058 

40.084 

148.13 

30 

30.027 

30.027 

18.410 

67.27 

35 

35.012 

35.012 

8.463 

30.70 

40 

40.006 

40.006 

3.996 

14.07 

Tropical: 

20 

20.137 

20.137 

95.154 

376.49 

25 

25.058 

25.058 

40.450 

1 50.46 

30 

30.026 

30.026 

18.306 

66.27 

35 

35.012 

35.012 

8.600 

30.37 

40 

40.006 

40.006 

4.181 

14.47 
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either  as  a  function  of  star  direction  and  satellite  position  only  or 
approximately  as  a  function  of  refraction  and  one  particular  component  of 
satellite  position  with  very  little  error. 

Figure  2-5  illustrates  how  hv  is  related  to  star  direction  and  satel¬ 
lite  position.  The  satellite  position  vector  (£)  can  be  expressed  as  the 
sum  of  two  orthogonal  vector  components  as  follows: 

r*fr.u)u  +  ( r  •  u  ]  u  (2-19) 

—  v—  —s'  — s  —up'  —up 

where  Ug  is  a  unit  vector  in  the  known  star  direction  and  Uup  is  a  unit 
vector  normal  to  Ug  and  in  the  plane  defined  by  Ug  and  £,  i.e. 

u  *  Unit  fu  x  (r  x  u  11  (2-20) 

—up  L— s  — s'J 

Note  in  Figure  2-5  that  hv  can  be  expressed  in  terms  of  the  magnitude  of 
the  second  vector  component  of  £  of  Equation  2-19  as  follows: 

i 

i 

h  =>  r  •  u  -  r  ( 2-21  ) 

v  —  —up  e 

Equation  2-21  therefore  gives  the  relationship  between  hv  and  satellite 
position  without  any  consideration  being  made  for  refraction. 

Figure  2-6  illustrates  how  hv  may  also  be  expressed  in  terms  of  the 
refraction  angle  (R),  the  apparent  ray  height  (ha),  and  the  magnitude  of 
the  other  orthogonal  component  of  satellite  position.  Using  triangle  BSE, 
the  vacuum  tangent  height  can  be  expressed  as  follows:  . 


h 


h  +  a  -  d  tan(R) 


(2-22) 


where  d  is  the  distance  of  the  satellite  from  the  cone  base  and  a  is  the 


small  distance  previously  indicated  to  be  negligible.  Note  that  d  is  equal 
to  the  magnitude  of  (£  •  u^ )  so  that  Equation  2-22  may  be  expressed  as 
follows: 

h  =  h  +  a  -  Ir  •  u  I  tan(R)  (2-23) 

v  a  | —  — s  I 

Equation  2-23  therefore  gives  the  relationship  between  hv  and  the  refrac¬ 
tion  cone.  In  addition,  this  equation  may  also  be  regarded  as  the 
measurement  equation  for  satellite  navigation  since  its  relates  hv  to  the 
measured  refraction  and  the  apparent  height  which  may  be  derived  from  the 
measured  refraction  for  a  known  atmospheric  density  profile.  The 
relationship  between  refraction  and  apparent  height  is  independent  of 
satellite  position  and  may  be  represented  by  the  following  general 
expression  which  expresses  ha  as  a  function  of  R  and  the  atmospheric 
density  profile  (denoted  merely  as  p ) : 

h  =  h  (R,  p)  (2-24) 

a  a 

Since  the  atmospheric  density  profile  is  assumed  to  be  known  in  the  case  of 
satellite  navigation,  the  apparent  height  for  a  given  refraction  can  be 
obtained  analytically  or  by  interpolating  a  table  of  precomputed  matching 
values  of  ha  and  R.  In  the  case  of  an  exponential  atmosphere,  the 
density  profile  is  exponential  by  definition  and  an  approximate  analytical 
expression  of  ha  (see  Appendix  A)  is  as  follows:  ' 


h  (R,  p  ) 
A  O 


h  -  H  ln(R)  +  H  In  1 

r2”  re_ 

r  H  r  -| 

kQ  )  p  - 

1  2 

|  +  R 

6 

°  1 

° L  h  J 

2  TT 

(2-25) 
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where  p0  is  the  density  at  some  arbitrary  reference  height  hQ,  and  H  is 


the  fixed  density  scale  height.  This  expression  will  be  used  later  to  in¬ 


dicate  the  effect  on  hv  of  errors  in  R  and  p •  For  the  moment,  however, 


the  function  in  Equation  2-24  is  assumed  to  be  that  for  any  known  arbitrary 


atmosphere,  and  when  substituted  for  h-  in  Equation  2-23  yields  the 


following: 


h  =  h  (R,  p )  +  a  -  r  •  u  tan(R) 
v  a  —  — s 


(2-26) 


This  expression  for  hv  provides  the  key  for  how  a  refraction  measurement 


be  used  to  update  satellite  position.  As  previously  shown  in  Equation 


2-19,  the  satellite  position  vector  can  be  expressed  as  the  sum  of  two 


orthogonal  vector  components  as  follows: 


r  *  fr  <  u  )  u  +  fr  •  u  lu 
—  —SJ—S  -up-'  —up 


(2-27) 


where  u_  is  the  unit  star  vector  and  u„n  is  a  unit  vector  normal  to 


Ug  and  in  the  plane  defined  by  Ug  and  r.  Although  u„p  is  normal  to 


the  star  direction,  and  not  the  cone's  surface,  the  difference  in  direction 


is  negligible  since  the  refraction  angles  are  small.  From  Equation  2-21 


the  magnitude  of  the  second  vector  component  of  Equation  2-27  can  be 


expressed  as  follows: 


r  •  u  =  h  +  r 
—  —up  v  e 


( 2-28 ) 


Substituting  the  expression  for  hv  from  Equation  2-26,  neglecting  the 


small  distance  a,  and  assuming  R  «  tan(R),  yields: 


r  *  =  h  (R,p)  -  r  •  u  R  +  r 

—  —up  a  —  — s  e 


( 2-29 : 
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which  when  substituted  into  Equation  2-27  yields  the  final  expression  re¬ 
lating  the  satellite  position  vector  to  the  refraction  angle  and  atmospher 
ic  density: 


r  =  f r  •  u  )  u  +  \ h  (R,p)  -  |r  •  u  I  R  +  r  1  u  (2-30) 

~  —s ;  -s  L  a  |—  — s  I  el  —up 

Here  it  is  seen  that  only  the  position  component  in  the  direction  of  unp 
would  be  updated  by  a’ refraction  measurement. 

2.4.4  Error  Relationships  Between  Satellite  Position,  Refraction  Angle  and 
Atmospheric  Density 

A  good  indication  of  the  effect  of  errors  in  refraction  measurement  and 
atmospheric  density  modeling  on  satellite  position  update  can  be  obtained 
by  considering  the  case  of  an  exponential  atmosphere  which  agrees  with  the 
1976  U.S .  Standard  Atmosphere  at  25  kilometers  altitude.  Substituting  the 
expression  for  ha(R,p)  from  Equation  2-25  into  Equation  2-29  and 
differentiating  gives: 


d(r 


u  1 

— up; 


f  H  r  1 

r  •  u  dR  + 

!  —  — s  1 

e 
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dR  ( 2-31  ) 


where  •  Ug |  is  the  distance  of  the  satellite  from  the  cone  base  which 
may  be  approximated  as: 


1  l 

12.  2  \ 

JL  •  “  ' 

(  r  -  (6400)  1 

(2-32) 


Inspection  of  Equation  2-31  indicates  a  relatively  simple  relationship 
between  the  satellite  position  error  (  d  (_r •  u^p ) ) ,  refraction  measurement 
error  (dR),  and  density  error  (dpQ).  Note  that  the  first  two  error  terms 
on  the  right  of  Equation  2-31  are  proportional  to  the  percentage  error  in 
density  or  refraction;  the  third  term  is  proportional  to  the  refraction 
error  and  the  distance  to  the  cone  base;  and  the  last  term,  which  happens 
to  be  very  small,  is  proportional  to  only  the  refraction  error. 

At  25  kilometers  altitude.  Tables  2-1  and  2-3  give  the  following  values 
for  the  U.S.  Standard  Atmosphere: 

H  =  6.366  km 
p0  =  40.084  gm/m^ 

R  =  148.1  arc  second  (0.000718  rad) 

For  a  one  percent  error  in  our  knowledge  of  density  at  25  kilometers, 
Equation  2-31  indicates  that  this  will  cause  a  satellite  position  error  of 
63.7  meters.  Likewise,  a  one  percent  error  (1.48  arc  seconds)  in 
refraction  measurement  for  a  satellite  at  an  altitude  of  1000  kilometers 
will  cause  an  overall  satellite  position  error  of  89.4  meters  (i.e.,  63.7 
meters  for  the  percentage  error  term  and  25.7  meters  for  the  last  two  terms 
of  Equation  2-31). 

2.4.5  Adoption  of  hv  as  Navigation  Measurement 

It  is  interesting  to  note  that  the  satellite  position  error  given  by 
Equation  2-31  can  also  be  interpreted  as  an  equivalent  error  in  the 


indicated  vacuum  tangent  height  since 


w  w  jtwijiuu*  a  i-i  uv  **  in  j*. 


h  =  r  •  u  -  r 
v  —  —up  e 


(from  Equation  2-28) 


dh  =  d  ( r  •  u  ") 
v  v—  -up-' 
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e 


0) 


dh 

v 
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H  r  -| 

__ _ e 

2  n  - 


(2-33) 


For  the  purposes  of  the  present  study.  Equation  2-33  suggests  that  it 
would  be  simpler  to  consider  hv,  rather  than  R,  as  the  parameter  being 
measured  since  it  conveniently  includes  the  error-effects  of  density  and 
refraction  measurement,  and  permits  one  to  express  the  overall  error  in 
common  units  of  meters.  This  approach  was  taken  in  the  present  study  and 
very  little  comment  is  made  in  the  later  sections  of  this  report  about  the 
relative  contributions  of  each  error  source  since  this  is  considered  to  be 
somewhat  arbitrary  in  a  general  study  of  this  type.  Most  of  the  navigation 
performance  data  generated  in  the  current  study  was  for  a  baseline  measure¬ 
ment  error  of  70  meters  in  hv,  although  some  data  was  generated  for  other 
values.  An  indication  of  the  error  in  hv  for  various  values  of  refrac¬ 
tion  measurement  error  combined  with  a  one  percent  error  in  density  is 
provided  in  Table  2-4  for  three  different  satellite  altitudes.  The  results 
were  generated  using  Equation  2-33.  Assuming  a  density  error  of  one  per¬ 
cent,  it  is  seen  that  the  baseline  value  of  70  meters  adopted  for  the  i 

current  study  is  representative  of  what  one  could  expect  with  a  refraction 
measurement  error  of  0.5  arc  second  at  a  satellite  altitude  of  1000  kilo¬ 
meters,  or  with  a  refraction  measurement  error  of  0.1  arc  second  at  geosta¬ 
tionary  altitude.  It  should  be  noted  that  the  star  trackers  on  the  satel-  .  < 

! 

f 

lite  HEAO-2  demonstrated  a  refraction  measurement  accuracy  of  better  than  j 

< 

0.5  arc  second.  ‘ 
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Table  2-4.  Equivalent  hv  Error  for  Various  Refraction  Measurement  Errors 

Combined  with  a  One  Percent  Density  Error.  (For  Observation  at 
25  km  Grazing  Altitude) 


Satellite 

Altitude 

(km) 

Density 

Error 

(%) 

Refraction 
Error 
(arcsec ) 

Equivalent  hv 

Error  (m)  For 

RSS 

po 

R 

(m) 

1000 

1 

1 .0 

64 

60 

88 

.5 

30 

71 

.1 

6 

64 

20,242 

1 .0 

168 

180 

(GPS) 

.5 

84 

106 

.1 

1  7 

66 

35,866 

1 .0 

276 

283 

(Geo) 

.5 

1  38 

152 

.1 

28 

70 
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2.4.6  Navigation  with  Dispersion  Measurements 

The  dispersion  represents  the  difference  in  refraction  between  two 
wavelengths  of  light.  According  to  Equation  2-6  the  refraction  at 
wavelength  X  for  an  exponential  atmosphere  is  approximately 


i 


R(A  )  *  k(A)  p  f  2nb  (r  +  h  ))  2 
g  '  e  g  * 


(2-34) 


where  pg  is  the  density  at  grazing  height  hg,  and  b  is  the  inverse 
density  scale  height  (i.e.  b  =  1/H).  Since  hg  <<  re,  the  refraction 
may  be  expressed  fairly  accurately  as 


R(A)  «  k(A)  p  f  2ifb  r  )  2 
g  v  e ' 


(2-35) 


Using  subscripts  b  and  r  to  denote  the  blue  and  red  wavelengths,  and 
denoting  p^  and  pr  as  the  densities  at  the  grazing  heights  for  the  blue 
and  red  rays,  the  dispersion  between  the  two  rays  is 


6  =  R.  -  R 
b  r 


«  =  (2*b  rj  2  (kfa  pb  -  kr  pr) 
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-b  (h  -  h  ) 

6  =  Rr  [K  e  gb  gr  -  1] 


(2-37) 


where  K  =  kb/kr. 

A  particular  value  of  5  sometimes  useful  in  the  present  instance  is 
that  value  commonly  referred  to  as  the  local  dispersion  angle  (6^).  This 
is  the  dispersion  between  the  red  and  blue  rays  for  the  same  grazing  height 
where 


=  Rr<  K-1 )  (2-38) 

It  is  important  to  note,  however,  that  a  satellite  cannot  observe  5^ 
since  the  red  and  blue  rays  for  the  same  grazing  height  will  diverge  after 
passing  through  the  atmosphere  and  are  sufficiently  separated  by  the  time 
they  reach  a  satellite  that  they  cannot  be  observed  simultaneously  with 
reasonable  size  optics.  Likewise,  a  satellite  cannot  observe  the 
dispersion  given  by  Equation  2-37  except  for  the  case  where  the  red  and 
blue  rays  converge  to  a  point  at  the  satellite  as  shown  in  Figure  2-7. 

Note  in  Figure  2-7  that  the  grazing  height  of  the  blue  ray  must  be  above 
that  of  the  red  ray  since  the  refraction  is  larger  for  the  blue  ray. 

To  guarantee  that  a  satellite  can  observe  the  dispersion  given  in 
Equation  2-37,  a  constraint  must  be  placed  on  the  difference  in  ray  grazing 
heights  shown  in  that  equation.  This  constraint  will  be  a  function  of  the 
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Fiqure  2-7.  Dispersion  geometry 


satellite  distance  from  the  horizon.  In  the  following  derivation,  the 


required  difference  in  grazing  heights  is  determined  by  first  expressing 
the  grazing  height  difference  in  terms  of  the  corresponding  difference  in 
apparent  heights.  Afterwards,  the  desired  constraint  is  placed  on  the 
apparent  heights  so  that  ray  convergence  occurs  at  the  satellite: 

According  to  Equation  2-18  the  grazing  height  difference  can  be 
expressed  in  terms  of  the  apparent  height  difference  as  follows: 


h  -  h 


h  -  k  p,r  h  -kpr 
b^be  ar  rKre 


1  +  kb  pb 


1  +  k  p 
r  K  r 


Since  the  denominator  terms  kb  pb  and  kr  pr  are  never  any  larger 
than  about  2  x  10-5  for  altitudes  above  20  kilometers,  they  can  be 
neglected  and 


Ah  -Ah  -  r  (k,  p,  -  k  p  ) 
g  a  ebb  r  r 


where 


Ah  =  h  -  h  and  Ah  =  h  -  h 
g  g,  g  a  a.  a 

b  r  b  r 


Substituting  for  p.  from  Equation  2-36 
b 


Since 


Ah  -  Ah  -k  p  r  (Ke^^g-i) 
g  a  r  r  e 


(2-39) 


e  ^^g  «  1  -  bAh  for  bAh  <<  1 

g  g 
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Equation  2-39  may  be  expressed  as 
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(2-41  ) 


To  determine  the  required  constraint  on  the  apparent  height  difference,  use 
is  first  made  of  Equation  2-15  to  obtain 

Aha  =  ha  "  ha  =  r3  [sin(Zb)  "  sin(Zr)]  (2_42) 

%  r 

where  Zb  and  Zr  are  the  zenith  angles  shown  in  Figure  2-7,  Note  in 
Figure  2-7  that  Equation  2-42  may  be  expressed  in  terms  of  the  angles  Eb 
and  Ej.  instead  of  the  zenith  angles  as  follows 

Ah  =  r  cos(E  ]  -  cosfE  ) 
a  s  v  h'  v  r' 


Aha  «  2  rs  sin  j(  Efa  ♦  Ej  sin^  -  E^) 


Since  <5  =  Ej.  -  Eh  and  E?  -  Eb 


Ah  •  6  r  sinf E  )  (2-43) 

a  s  v  r* 


Note  that  Equation  2-43  could  have  been  obtained  by  inspection  of  Figure 
2-7  since  Aha  is  the  distance  subtended  by  6  at  a  satellite-to- 
horizon  distance  of  rg  sin(Er).  Substituting  Equation  2-43  into 
Equation  2-41  yields  the  desired  constraint  on  the  grazing  height 
difference 


Ah 


5  r  sin(E  )  -  k  p  r  ( K— 1  ) 
s _ r _ r  r  e _ 

1  -  K  b  k  p  r 
r  Kr  e 


(2-44) 


Using  the  approximation  of  Equation  2-40  in  Equation  2-37,  substituting  the 
expression  for  Ahg  from  Equation  2-44  into  Equation  2-37,  and 
substituting  the  equivalent  expression  for  pr  in  accordance  with  Equation 
2-35,  yields  the  following  desired  expression  for  the  dispersion  observed 
by  a  satellite: 


6  -  6*/C 

where  6^  is  given  in  Equation  2-38  and 


(2-45) 


_1_ 

C  =  1  +  K  b  [r  sin(Er)  -  (r  /2wb)  2  ]  (2-46) 

An  indication  of  how  well  Equation  2-45  indicates  the  dispersion  for  a 
realistic  atmosphere,  such  as  the  1976  U.S.  Standard  Atmosphere,  is  shown 


in  Table  2-5  where  the  results  are  shown  along  with  those  generated  by 


Table  2.5.  Observed  Dispersion  Between  0.35  and  0.7  Microns  for  U.S.  Standard  Atmosphere  and  Two 

Satellite  Altitudes 
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accurate  ray  tracing^11).  The  values  of  b  and  Rr  used  to  generate  the 
results  with  Equation  2-45  are  those  given  for  the  1976  U.S.  Standard 
Atmosphere • 

Note  that,  for  a  given  grazing  height,  the  magnitude  of  the  factor  C 
increases  with  satellite  altitude,  thereby  reducing  the  magnitude  of  the 
observed  dispersion.  However,  as  previously  noted,  the  dispersion  can  be 
measured  very  accurately  with  a  two-color  ref ractometer .  The  potential 
accuracy  of  this  device  is  a  function  of  various  parameters  such  as 
telescope  aperture  size,  star  brightness,  and  optical  transmission 
characteristics.  Tests  conducted  by  the  University  of  Maryland  have 
indicated  an  accuracy  better  than  0.001  arc  second. 

As  previously  indicated  in  Section  2.2,  the  dispersion  measurement  can 

be  converted  to  an  equivalent  refraction  measurement.  Equation  2-45 

provides  a  fairly  good  indication  of  what  the  relationship  would  be  between 

the  dispersion  for  two  wavelengths  (red  and  blue)  and  the  refraction  Rr 

for  the  red  wavelength.  It  will  be  recalled  that  the  locus  of  constant 

refraction  may  be  visualized  as  an  earth-based  cone  extending  out  into 

space,  thereby  permitting  a  relatively  simple  relationship  between 

refraction  and  satellite  position.  Such  is  not  the  case  for  dispersion 

since  the  locus  of  constant  dispersion  is  a  more  complicated  surface  in 

space.  However,  by  converting  the  dispersion  measurement  to  an  equivalent 

refraction  measurement,  one  can  make  use  of  the  latter's  relatively  simple 

cone  geometry  to  relate  the  dispersion  measurement  to  satellite  position. 

« 

A  fairly  good  indication  of  the  equivalent  refraction  is  obtained  by 


inverting  Equation  2-45: 


5 


(2-47) 


\  —  i 

r  K-1  -  KbS  [  rg  sin(Er)  -  (re/2*b)  2] 


The  error  relationship  between  dispersion  and  its  equivalent  refraction 


is  obtained  by  differentiating  either  Equation  2-47  or  2-45: 


dRr  -  as  -  jrr  as 


(2-48) 


where  G  is  the  denominator  of  Equation  2-47  and  C  is  the  denominator  of 


Equation  2-45. 


Substituting  the  last  expression  in  Equation  2-48  for  dR  in  Equation 


2-33  yields  the  following  relationship  between  the  errors  in  dispersion, 


atmospheric  density  modelling  and  vacuum  tangent  height: 


■  *r  *  [("  V2*)  2  -  H/Rr  -  |r  •  2.1]  ib 


d5  (2-49) 


Substituting  b  =  1/H  and  assuming  rssin(Er)  ■  r_  •  Ug  ,  one  can  use 


the  relationship  of  Equation  2-46  to  re-express  the  bracketed  expression  of 


Equation  2-49  so  that 


dh  = 
v 


1  r i  .  c  (c  +  K-1) 

■  b  p  dpo  K  R  (K  -  1  )  d6 
-  o  r 


(2-50) 
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The  sensitivity  of  dispersion  measurement  to  satellite  position  (e.g., 


hv)  is: 


K  b  R  (K  -  1 ) 
_ r _ 

C2  (C  +  K  -  1 ) 


(2-51  ) 


Figure  2-8  shows  how  well  the  sensitivity  values  generated  by  Equation  2-51 
compare  with  those  obtained  by  differentiating  the  ray  trace  data  with  a 
"5-point  derivative"  formula  for  a  satellite  at  two  different  altitudes. 
Note  that  the  ray  trace  curves  in  Figure  2-8  contain  some  high  frequency 
oscillation  due  to  differentiation  of  the  discrete  data  associated  with  a 
multi-layered  atmospheric  model  consisting  of  72  levels  between  14  km  and 
100  km  altitude.  Also  note  that  the  maximum  sensitivity  for  the 
higher -altitude  satellite  occurs  at  a  higher  grazing  height  than  that  for 
the  lower -altitude  satellite.  By  differentiating  Equation  2-51  it  is  found 
that  the  maximum  sensitivity  occurs  when 

C2  +  (2--|)C  +  K-1  =  0  (2-52) 


For  the  wavelengths  being  considered  in  this  report  (X  =  0.35  and  0.7 
microns),  Equation  2-52  is  satisfied  when 

C  =  -0.02485,  1.5061  t  (2-53) 

where  only  the  latter  value  is  realizable.  Substituting  this  value  of  C 
into  Equation  2-46  yields  the  associated  value  of  red -wave length  refraction 
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in  radians  (denoted  as  R  ) : 
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Figure  2-8.  Dispersion  angle  sensitivity  vs  grazing  ray  tangent  height 
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(2-54) 


=  0.48787/t  fr s  sinfE^)  -  (r^wb)2 


If  an  exponential  atmosphere  is  assumed,  the  grazing  height  (hg),  where 
maximum  dispersion  sensitivity  occurs,  may  be  determined  by  substituting 
Rr  for  R  in  Equation  A-6: 


where  pQ  is  the  density  specified  at  some  arbitrary  height  (hQ)  in  the 
exponential  atmosphere.  Using  pQ  =  40.084  g/m^  at  hQ  =  25  km  from 

I 

the  1976  U.S .  Standard  Atmosphere,  the  values  of  hg  computed  for 
satellite  altitudes  of  1,000  and  20,000  km  are  23.8  and  36.3  km. 
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SECTION  3 


ERROR  ANALYSIS 

3.1  Introduction 

3.1.1  Objectives 

In  choosing  an  error  model,  the  following  objectives  were  considered. 

1)  Determination  of  the  navigator's  steady  state  performance 
(position  error  standard  deviations)  assuming  certain  baseline 
parameters  for  its  measurement  capabilities. 

2)  Sensitivity  analysis  of  the  navigator  to  changes  in  the  number 
and  direction  of  star  sightings  and  to  changes  in  atmospheric 
modeling. 

3)  Estimation  of  the  navigator's  settling  time  (time  for  the 
position  error  standard  deviations  to  approach  their  nominal 
steady  state  values). 

Review  of  these  goals  led  to  the  selection  of  a  linearized  Kalman  filter. 
This  filter  is  fundamentally  the  same  algorithm  as  would  be  incorporated  i 
the  estimator  of  the  actual  navigation  system. 

3.1.2  Filter 

The  Kalman  filter  is  a  numerical  algorithm  which  recursively  estimates 
the  current  state  of  a  system,  given  external  measurements,  z.  The 
algorithm  incorporates  models  of  the  system  and  measurement  dynamics,  and 
assumed  statistics  of  system  noises  and  measurement  errors,  to  produce 
estimates  which  minimize  the  mean-squared  errors  (variances)  of  the  state 
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variables.  If  the  system  is  linear  and  its  errors  normally  distributed, 
then  the  Kalman  filter  is  truly  the  optimal  estimator  (see  Gelb^14^). 

The  discrete  linearized  Kalman  filter,  the  estimator  chosen  for  this 
analysis,  is  shown  in  Table  3-1.  The  non-linear  system  of  orbital  motion 
is  linearized  about  a  reference  trajectory  which  nominally  describes  the 
time  history  of  the  state  variables;  the  state  vector  and  covariance  matrix 
are  then  propagated  along  this  path;  this  filter  provides  a  nearly  optimal 
estimate  of  the  non-linear  system  with  considerably  less  computation  than 
more  accurate  non-linear  filters. 

A  fundamental  component  of  the  Kalman  filter  is  the  covariance  matrix, 
P^.  The  diagonal  terms  of  Pi  are  the  error  variances  of  the  state 
variables;  off-diagonal  terms  are  the  covariances  between  pairs  of  state 
variables.  The  elements  of  the  covariance  matrix,  especially  the  diagonal 
terms,  provide  a  yardstick  for  determining  how  precisely  the  filter  is 
estimating  the  state. 

A  very  valuable  property  of  the  Kalman  filter  is  that  the  error 
covariance  matrix  can  be  propagated  and  updated  without  estimating  the 
state  itself.  This  represents  a  significant  saving  in  computer  time 
compared  to  other  techniques.  Such  covariance  analyses  have  become  widely 
accepted  tools  in  analysing  system  performance. 

3.2  System  Description 

3.2.1  States 

The  linearized  state  of  the  navigator  may  be  described  b^  seven 
elements:  three  components  of  position  error,  three  components  of  velocity 

error  and  a  scalar  measurement  bias  error.  These  elements  form  the  state 
vector 
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with  the  error  covariance  matrix 

P  =  E{x  xT}  (3-2) 

where  6r_  and  6v^  define  the  errors  in  the  trajectory  estimated  by  the 
navigator.  Their  units  are  kilometers  and  kilometers/second,,  respectively. 

The  scalar  measurement  bias  error  was  introduced  as  am  additional 
variable  in  the  above  state  to  account  to  some  extent  for  bias-type  errors 
in  atmospheric  modeling  and  refraction  or  dispersion  measurements.  The 
addition  of  this  variable  to  the  estimated  state  permits  the  navigator  to 
remove  its  effect  from  the  estimated  position  and  velocity.  It  should  be 
noted,  however,  that  the  addition  of  this  state  variable  (and  its 
associated  variance)  adds  complexity  and  randomness  to  the  system  and  so 
increases  the  size  of  the  position  and  velocity  error  variances  compared  to 
the  case  where  no  bias  error  is  present.  In  the  present  study  this  bias 
error  is  treated  as  the  overall  bias  error  in  the  measurement  of  vacuum 
tangent  height  resulting  from  systematic  errors  in  atmospheric  modeling  and 
refraction  or  dispersion  measurements.  In  Section  2.4.5  the  relationship 
between  errors  in  vacuum  tangent  height,  atmospheric  modeling  and 
refraction  measurement  are  given  for  an  exponential  atmosphere.  This 
relationship  indicates  that  a  more  elaborate  representation  of  systematic 
error  could  have  been  included  in  the  present  analysis  (i.e.  separate 
systematic  errors  for  density,  refraction  and  density  scale  height). 
However,  for  the  purposes  of  this  preliminary  investigation,  it  was  felt 
that  the  use  of  a  bias  error  in  vacuum  tangent  height  was  sufficient. 


3.2.2  Coordinate  Systems 

The  satellite's  position  and  velocity  may  be  described  in 
earth-centered  inertial  (ECI)  coordinates.  The  z  axis  points  north  along 
the  Earth's  polar  axis,  the  x  axis  points  in  the  direction  of  the  vernal 
equinox  and  the  y  axis  completes  a  right-handed  coordinate  system.  This 
system  is  the  one  ultimatly  used  to  propagate  and  update  the  covariance 
matrix. 

While  the  ECI  coordinate  system  is  convenient  for  computation,  the 
position  uncertainties  and  update  directions  are  more  readily  visualized  in 
terms  of  the  satellite-centered  coordinate  system  TNR  (Tangential,  Normal, 
Radial)  shown  in  Figure  3-1.  In  this  system  u^  points  toward  the  earth 
along  the  local  vertical,  ujj  is  normal  to  the  orbital  plane  and  points  in 
the  opposite  direction  of  the  satellite's  orbital  angular  momentum  vector; 
u<r  is  orthogonal  to  uj*  and  ujj  and  is  directed  in  the  same  sense  as 
the  satellite's  velocity  vector.  In  Figure  3-1  the  direction  of  a  star 
(Ug )  undergoing  refraction  can  be  specified  in  TNR  coordinates  by  an 
azimuth  angle  (p)  and  elevation  angle  (0  )  where  the  latter  is  also  shown  in 
Figure  3-2. 

The  direction  of  satellite  position  update  for  the  star  in  Figures  3-1 
and  3-2  is  given  by  the  unit  vector  Uyp  where 


u  =  cos  0  f cos  A  u,  +  sin  A  u  1  -  sin  0  u 
-up  ^  y  -T  Y  — N'  —i 


(3-3) 


The  relationship  between  the  ECI  and  TNR  coordinate  systems  is  given  by 


TNR  TNR  ECI 

r  =  c  T  r 
—  ECI  — 


(3-4) 
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The  transformation  for  the  full  state  vector  is  given  by 


TOR  *  EC I 
x  =  C  x 


(3-6) 


where 
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Note  that  the  measurement  bias,  a  scalar,  is  not  affected  by  the  trans¬ 
formation  of  coordinates. 


3.3  Propagation 

3.3.1  State  Transition  Matrix 

The  state  transition  matrix  $(t,tQ)  permits  the  direct  and  simple 
linear  calculation  of  the  state  vector  at  time  t  given  the  state  at  time 
to,  provided  that  the  system  is  unforced  during  that  time  interval. 


Thus, 


jc(t)  =  * (t,tQ)x(t0) 


(3-8) 


Similarly,  the  error  covariance  matrix  P  at  time  tg  may  be  propagated  to 
time  t  using 


P(t)  =  #(t,t0)P{tQ)*  (t,tQ) 


(3-9) 


Let  x  be  the  navigator's  reduced  state  vector  (position  and  velocity) 
and  $*(t,tg)  be  the  corresponding  transition  matrix  between  time  t  and 
tg •  Let 


G(t) 


r=r  _ 
- ref 


(3-10) 


where 


£ 

£ref 


local  gravitational  acceleration 
reference  trajectory 


and 


F(t) 


(3-11  ) 


Then  $*(t,tg)  is  the  solution  of  the  matrix  differential  equation 


d_ 

dt 


♦  *(t. 
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F(t)**(t,  tQ) 


(3-12) 


subject  to  the  initial  conditions 
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These  equations  may  be  solved  in  closed  form,  using  methods  described  by 


Goodyear ( 1 . 


The  transition  matrix  for  the  full  state  vector  x  is  related  to  the 


reduced  state  transition  matrix  by 


*  i 

*  (t,t.)  ,  0 


*<t.t0)  = 


0 '  I  U6x1 
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°6x  i  ;  1 


(3-14) 


Note  that  the  estimate  of  the  measurement  error  bias  remains  constant 


during  the  transition. 


3.3.2  Process  Noise 


Disturbing  accelerations  strongly  affect  the  long  term  behavior  of  an 


orbit.  Such  disturbances ^ 16)  are  introduced  by  the  higher-order  gravity 


field  terms,  atmospheric  drag,  solar  pressure,  etc.  To  correctly  model 


these  effects  in  the  present  study  would  be  computationally  expensive  and 


somewhat  unnecessary.  Therefore,  these  disturbances  were  treated  as 


process  noise  driving  the  unforced  continuous  system,  i.e. 


X  =  Fx  +  w 


(3-15) 


where 


x  =  system  state  vector 


w  =  unmodeled  disturbance  vector 


-  64  - 


since  the  noises  (disturbing  accelerations)  directly  act  on  the  velocities 
only.  These  noises  are  treated  as  white  and  unbiased,  with  the  spectral 
density  Q  for  the  discrete  case  given  as 


T 

E{w  w  } 


(3-17) 


The  propagation  of  the  error  covariance  matrix  is  similarly  disturbed 
by  an  additive  term 


P.  ,  =  *(t.  ,,t.  )P.$  (t.  ,,t.  )  +  Q. 
i+1  i+I  i  i  1+1  l  l 


(3-18) 


where  Qj^  is  the  covariance  matrix  of  the  process  noise. 

Formally,  the  transition  from  the  ncise  spectral  density  matrix  to  the 
noise  covariance  matrix  is  given  by 

=  /ti  +  1*  (ti+1  ,T  )Q(t  )$T(ti+1  ,T  )dT  (3-19) 

i 


However,  for  small  time  step  At  this  equation  can  be  approximated  by 

% 

Q  =  QAt  (3-20) 


(see  Schweppe ( 1 7 ) ) . 


3.4  Measurement  Incorporation 

The  linear  Kalman  filter  requires  that  the  system  measurements  z_  be 
linear  combinations  of  the  current  state  x  and  a  random  measurement  error, 
e^ 

z .  =  H.  x.  +  e .  ( 3-21  ) 

—l  l—i  ~ l 

The  error  covariance  matrix  P^  is  updated  using  only  the  measurement 
matrix  H^. 

P.  (  +  )  =  f  I  -  K.  H.  Ip.  (-)f  I  -  K.  H.  1T  +  K.  R.  KT  (3-22) 
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where 


(3-23) 


and 


R. 

l 


(3-24) 


What  follows  are  derivations  of  the  matrices  required  for  the 
covariance  update. 

3.4.1  Measurement  Matrix 

By  observing  starlight  refraction  through  the  atmosphere,  ’the  navigator 
effectively  measures  the  vacuum  tangent  height  with  an  unknown  measurement 
bias  b.  Thus 

z  =  h  +  b  ( 3-25 ) 

v 
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The  bias  and  its  variance  may  not  be  the  same  everywhere  and  some 


consideration  is  given  to  its  treatment  in  Section  4. 3. 1.3. 
According  to  Equation  2-21 


h  =  r  •  u  -  r 
v  —  — p  e 


where 


r  =  satellite  position 


re  =  earth  radius 


and  u„p  is  the  update  direction  (Equation  2-20)  given  by 


u  =  Unit  [  (u  x  (r  x  u  ] 
—up  1  — s  —  — sJ 


By  definition 


where  x  is  the  non-linearized  state  vector 


x  =  v 


(3-26) 


(3-27) 
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Since  z  is  not  a  function  of  velocity 
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(3-30) 


Solving  for  the  measurement  gradient  with  respect  to  position: 
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(3-31  ) 


It  can  be  shown,  however. 


that  the  gradient  vector 


is 


normal  to  r,  i.e. 
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Therefore, 
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(3-33) 


3.4.2  Measurement  Error  Matrix 

Since  the  navigator's  measurement  is  a  scalar,  the  measurement  error 
matrix  R^  (see  Equation  3-24)  contains  only  the  single  variance 
associated  with  the  error  in  estimating  the  vacuum  tangent  height  hv, 

i.e. 

2 

=  <Ju 


R. 


(3-34) 


SECTION  4 


IMPLEMENTATION 


4.1  Programs 

Two  different  computer  programs,  arbitrarily  designated  as  SHAD1  and 
SHAD2,  were  developed  to  implement  the  Kalman  filter  described  in  Section 
3;  each  used  a  different  means  to  simulate  star  sightings  and  orbital 
trajectories.  By  comparing  the  results  of  SHAD1  with  those  of  SHAD2,  one 
can  judge  the  validity  of  the  assumptions  made  in  each. 

4.1.1  SHAD1 

SHAD1  used  the  the  following  approximations  in  its  implementation  of 
the  navigation  filter. 

1)  The  satellite's  orbit  is  Keplerian. 

2)  Times  and  directions  (in  satellite-centered  coordinates)  of 
measurement  stars  (those  appearing  close  to  the  horizon)  were 
randomly  generated. 

3)  The  covariance  matrix  was  propagated  and  updated  only  at 
uniformly  spaced  time  intervals. 

i 

SHAD1  did  not  utilize  a  catalog  of  star  positions  in  terms  of  celestial 
coordinates.  Rather,  the  program  selected  a  number  of  random  sighting 


times  where  the  satellite  would  observe  a  grazing  star.  For  each  sighting 


time,  an  independently  generated  random  star  direction  was  assigned  which 


placed  the  star  at  a  specified  tangent  height  when  viewed  from  the  satel¬ 
lite.  The  sighting  times  and  star  directions  were  generated  in  a  manner  which 
is  consistent  with  the  probability  distribution  expected  of  measurements 
taken  from  a  uniformly  distributed  star  field. 

4.1 .2  SHAD 2 

SHAD2  used  more  realistic  assumptions  for  the  navigator  than  SHAD1 . 

1 )  Any  trajectory  can  be  used,  provided  the  positions  and  velocities 
are  given  for  the  times  of  interest. 

2)  Star  directions  for  measurements  were  obtained  from  the  Yale 
University  Catalog  of  Bright  Stars. 

3)  The  covariance  matrix  was  propagated  between  irregularly  spaced 
time  intervals  corresponding  to  the  times  between  actual 
measurements. 

4.2  Progagation 

4.2.1  Trajectories 

Both  SHAD1  and  SHAD2  required  a  nominal  trajectory  (time,  position,  and 
velocity  of  the  satellite  at  frequent  intervals)  about  which  ‘to  propagate  the 
err.  r  covariance  matrix.  Auxiliary  programs  were  used  to  generate  the 
trajectory  files  for  the  two  filter  programs.  For  SHAD2,  use  was  made  of  the 
Goddard  Trajectory  Determination  System  (GTDS)^1®)  to  generate  satellite 


trajectories  perturbed  by  the  dominant  higher-order  gravitational  term  J2* 


4.2.2  State  Transition  Matrices 


Propagation  of  the  covariance  matrix  between  updates  was  accomplished 
with  state  transition  matrices  computed  for  the  time  steps  required.  A 
technique  developed  by  Goodyear ( 1 5 )  was  used  to  calculate  the  state 
transition  matrices  for  the  Keplerian  orbits  of  SHAD1 . 

Since  the  direct  calculation  of  state  transition  matrices  for  the 
perturbed  orbits  of  SHAD2  was  considered  to  be  very  costly,  a  method  of 
patched  conics  was  utilized  to  obtain  an  approximate  solution.  Insofar  as 
propagation  of  the  covariance  matrix  is  relatively  insensitive  to  small 
changes  in  the  state  transition  matrix,  this  approximation  was  considered 
to  be  satisfactory.  Given  a  description  of  a  perturbed  orbit  with  closely 
spaced  time  steps,  the  method  of  patched  conics  was  used  as  follows.  The 
instantaneous  position  and  velocity  of  a  satellite  describes  a  Keplerian 
orbit  which  is  tangent  to  the  satellite’s  actual  path  at  that  point.  Over 
a  short  time  interval  this  instantaneous  conical  path,  or  osculating  orbit, 
remains  very  close  to  the  perturbed  trajectory.  By  calculating  the 
osculating  orbits  for  all  known  points  in  the  perturbed  orbit,  one  can 
approximate  the  true  trajectory  with  one  consisting  of  a  series  of  conical 
orbits  patched  together ^ 1 9 ^ .  The  state  transition  matrices  along  these 
osculating  orbits  are  readily  calculated  with  the  technique  employed  in 
SHAD1 . 


4.2.3  Estimation  of  Process  Noise 

As  previously  indicated  in  Section  3.3.2,  a  process  noise'  was  included 
in  the  linearized  filter  to  statistically  account  for  the  un-modeled 
disturbing  accelerations  due  to  air  drag,  solar  pressure  and  the  higher 
order  gravitational  terms.  This  noise  was  estimated  by  first  sampling  the 


disturbing  accelerations  indicated  by  the  GTDS  program  at  discrete  times 


throughout  the  satellite  orbit  and  estimating  the  variances  of  the  samples 
in  TNR  coordinates.  A  spectral  density  matrix  Q  was  then  formed,  with  all 
of  its  non-zero  diagonal  elements  being  assigned  the  value  of  the  largest 
variance  component.  This  assured  that  the  probability  distribution  of  the 
process  noise  would  be  spherical  and  large  enough  to  conservatively  model 
the  system.  From  the  spectral  density  matrix,  the  process  noise  covariance 
matrix  was  calculated  as  shown  in  Equation  3-20. 

It  should  be  noted  that  the  process  noise  used  in  this  study  included 
the  effects  of  all  higher  order  gravitational  terms  of  the  Goddard  Earth 
Model  Number  9  (GEM9 )  except  for  the  dominant  oblateness  term  J2  whose 
effect  was  already  included  in  the  generation  of  satellite  trajectories  for 
SHAD2.  This  process  noise  was  used  in  the  generation  of  performance 
results  for  both  SHAD1  and  SHAD 2 .  Consequently,  the  performance  results  of 
SHAD1  are  those  for  the  case  where  the  precession  effect  of  J2  is  assumed 
to  be  non-existent,  while  those  of  SHAD2  do  include  its  effect.  An  actual 
navigator  would  probably  incorporate  the  dynamics  of  J2,  and  possibly 
some  additional  terms,  in  the  state  transition  matrix  and  treat  all  other 
disturbances  as  process  noise.  However,  for  the  purposes  of  the  present 
covariance  analysis,  it  was  felt  that  the  performance  results  would  be 
fairly  indicative  of  what  one  can  expect. 

4.3  Measurements 


4.3.1  Measurement  Parameters 
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4. 3. 1.1  Number  of  Stars 


The  number  of  star  sightings,  and  hence  the  number  of  covariance 
updates,  performed  during  an  orbital  revolution  is  primarily  dependent  on 
the  brightness  limits  of  the  navigator's  star  sensor.  As  the  navigator's 
sensitivity  is  improved,  more  stars  can  be  acquired  for  navigation.  This 
can  be  simulated  in  the  program  by  varying  the  number  of  star  positions 
made  available  for  covariance  updates. 

SHAD1  requires  the  user  to  enter  the  number  of  stars  in  the  sky  which 
are  detectable  (i.e.  of  sufficient  brightness)  to  the  navigator.  SHAD1  then 
estimates  the  average  number  of  stars  that  the  navigator  can  expect  to  view 
during  a  single  orbital  rotation  and  generates  the  appropriate  number  of 
star  sighting  directions. 

SHAD 2  employs  the  Yale  University  Catalog  of  Bright  Stars  to  describe 
its  star  field  and  uses  brightness  limits  to  decide  whether  or  not  an 
individual  grazing  star  will  be  sighted.  If  not,  the  program  will  continue 
covariance  propagation  through  the  time  that  a  star  sighting  might  have 
been  made . 

4. 3. 1.2  Azimuth  Span 

The  azimuth  span  specifies  how  much  and  what  part  of  the  horizon  will 

be  scanned  for  grazing  stars  (see  Figure  4-1).  Increasing  the  width  of 

this  span  places  more  stars  within  the  navigator's  field  of  view.  For 

applications  requiring  extreme  accuracy  this  span  should  encompass  360 

* 

degrees.  In  most  cases  the  azimuth  span  may  exclude  the  forward  section  of 
the  horizon  since  rising  stars  are  more  difficult  to  locate  before  they 
leave  the  atmosphere.  Stars  near  ±90  degrees  azimuth  will  take  much  longer 


i  ■*.  WAARiumji  rj  MKjrwrwi*jnuni  inur*  v*ywrj«  '^wwiwwjvuvwvwtww. 


to  go  into  refraction  than  those  close  to  the  orbital  plane.  In  the  time 
taken  to  observe  one  star  near  ±90  degress  it  might  be  possible  to  observe 
two  or  more  stars  closer  to  the  orbital  plane.  For  this  reason,  most 
simulations  use  am  azimuth  spam  between  135  and  225  degrees. 
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where  *  earth's  polar  axis 


and 


u  =  Unit  fu  x  (r  x  u  )1  (4-2) 
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Both  SHAD1  and  SHAD 2  use  the  grazing  latitude  to  model  the  measurements  of 
atmospheric  refraction  in  one  of  three  different  ways: 


1  )  Bias  Everywhere 

The  program  treats  every  measurement  as  if  it  contained  a  bias  which 
can  be  estimated.  Measurements  may  be  taken  over  any  part  of  the 
earth.  The  measurement  matrix  H  is  always  given  by 


r  T  i  < 
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(4-3) 


2)  No  Bias,  Measurements  Limited 

Refraction  measurements  are  limited  to  those  for  which  the  grazing 
latitude  falls  within  a  band  about  the  equator  given  by  ±\|»max  (see 
Figure  4-2).  Within  this  band  atmospheric  behavior  is  understood 
well  enough  to  model  measurements  without  a  bias.  Any  star  found  to 
graze  outside  of  this  band  will  be  ignored  for  measurements. 

Since  a  measurement  bias  is  not  observed  during  a  measurement,  the 


measurement  matrix  becomes 


Without  the  observation  of  a  bias,  no  cross  correlation  can  be  built 
up  in  the  covariance  matrix  between  the  measurement  bias  variance 
and  other  state  variances.  The  value  chosen  for  the  initial 
measurement  bias  variance  can  be  arbitrary  since  it  will  not  be 
changed  by  a  measurement  nor  will  it  affect  the  other  state 
variances.  Effectively,  the  state  of  the  system  is  reduced  by  one 
state  variable. 


Figure  4-2.  Measurement  restrictions  due  to  grazing  latitude 
In  the  cases  tested,  the  maximum  grazing  latitude  was  set  at  30 


degrees  for  this  option 


Refraction  measurements  are  made  everywhere  but  measurements  made 


outside  of  a  specified  equatorial  band  ¥max  are  treated  as 
containing  a  bias  while  those  inside  are  not.  For  an  individual  star 
sighting,  if  its  grazing  latitude  is  within  the  equatorial  band  then 
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otherwise 
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4.3.2  Simulated  Star  Field  (SHAD1 ) 


4. 3. 2.1  Description 

The  brighter  stars  visible  from  the  Earth  are  distributed  almost 
uniformly  throughout  the  sky.  This  makes  it  possible  to  greatly  simplify  the 
simulation  by  updating  the  filter  with  artificially  generated  star 
sightings.  The  artificial  star  field  should  yield  approximately  the  same 
navigator  performance  as  an  actual  star  field. 

These  random,  artificially  generated  star  positions  need  not  be  specified 
in  inertial  space.  Suppose  we  have  arbitrarily  assigned  a  star.  This  star 
must  be  placed  in  the  sky  so  that  its  light  will  graze  at  the  desired  tangent 
altitude.  For  convenience,  let  the  star  direction  be  given  in  TNR 
coordinates.  If  we  assume  a  spherical  earth,  the  horizon  elevation  angle  0 


is  specified  by  the  satellite  and  vacuum  tangent  altitude  (see  Figure  3-2):  a 
star's  position  can  thus  be  defined  given  a  sighting  azimuth 


An  artificial  star  field  can  be  created  by  randomly  assigning  azimuths 
and  a  schedule  of  sighting  times  (tantamount  to  specifying  positions  in  the 
known  reference  orbit).  If  desired,  this  data  can  be  converted  to 
represent  star  positions  in  ECI  coordinates. 

4. 3. 2. 2  Probability  Distribution 

To  simulate  a  realistic  star  field,  the  schedule  of  star  sightings 
(sighting  times  and  azimuths)  must  be  consistent  with  the  presumption  that 
bright  stars  are  uniformly  distributed  over  the  sky.  Therefore,  the 
probability  distributions  used  to  randomly  generate  sighting  times  and 
azimuths  must  be  such  that  they  conform  to  a  uniform  star  field  when 
converted  to  ECI  coordinates. 

m 

Circling  about  the  earth,  the  navigator  observes  stars  near  the  horizon 
moving  into  and  out  of  refraction.  If  its  orbit  is  circular  the  navigator 
has  a  constant  expectation  of  viewing  a  grazing  star  anywhere  in  the 
orbit.  In  this  case,  the  orbital  angular  positions  of  the  satellite  when 
star  sightings  occur  will  be  uniformly  distributed  between  0  and  360 
degrees.  These  angular  positions  can  be  readily  converted  into  times  when 
the  navigator  observes  a  star. 

As  the  navigator  scans  the  horizon  within  its  azimuth  span,  it  will,  in 
the  course  of  one  orbital  period,  survey  a  sector  of  the  celestial  sphere 
(see  Figure  4-3).  It  can  be  readily  shown  that  the  fraction  of  the 
celestial  sphere  scanned  is  given  by 

Fr  =  9  sin  $2  -  sin  ^  (4-7) 
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Figure  4-3.  Sector  of  the  celestial  sphere  scanned  by  navigator. 

where  $1  and  $2  are  l°wer  and  upper  limits  of  the  azimuth  span. 

Clearly,  if  we  specify  that  there  shall  be  N  stars  on  the  celestial  sphere 
of  sufficient  brightness  to  be  used  by  the  navigator,  then,  on  the  average 
only  N*  of  these  will  actually  be  usable  during  one  orbital  period  where 

* 

N  = 


Fr  N 


(4-8) 


The  distribution  of  star  sightings  as  a  function  of  azimuth  may  be 
derived  by  observing  that  the  probability  of  observing  a  star  within  a 
narrow  annular  segment  of  width  d$  is  proportional  to  the  area  of  that 
segment;  the  differential  area  depends  upon  its  azimuth.  It  can  be  shown 
that  the  probability  density  of  star  sightings  as  a  function  of  azimuth  is 


P(<|>  ) 


_ COS  it _ 

sin  $2  ~  sin 


(4-9) 


The  user  specifies  the  number  of  stars,  N,  which  he  wishes  to  place  in 
the  artificial  sky  and  the  azimuth  span,  and  $2*  of  the  sensor.  For 

a  given  orbit,  SHAD1  estimates  the  number  of  stars  that  the  sensor  can 
expect  to  find  within  its  azimuth  span  and  generates  that  number  of  star 
sightings  and  directions  for  each  orbit. 

4. 3. 2.3  Sighting  Times 

Sighting  times  were  created  by  generating  random  numbers  uniformly 
distributed  between  0  and  1 .  The  random  numbers  are  used  to  specify  the 
fraction  of  the  orbit  traveled  from  the  satellite's  initial  position.  For 
a  circular  orbit,  this  fraction  readily  gives  the  time  of  flight  since  the 
satellite  last  passed  its  initial  position.  Sighting  times  for  additional 
revolutions  were  obtained  by  adding  the  appropriate  number  of  orbital 
periods . 

4. 3. 2.4  Azimuth  Angles 

To  generate  random  azimuth  angles  such  that  they  have  the  probability 
distribution  given  in  Equation  (4-9),  use  was  made  of  the  probability 
function  Pr  ($  £$).  Let  us  consider  the  random  variable 
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s  =  Pr  {(f> 


(4-10) 


where  $  is  one  of  a  random  ensemble  of  azimuths  having  the  desired 
distribution.  By  definition,  s  is  uniformly  distributed  between  0  and  1. 
Furthermore,  every  value  of  s  uniquely  specifies  a  value  for  <j>  by 


=  sin  [s(<j>)  (sin  4> 2  -  sin  1 )  +  sin  <j> 


!♦!  <  t 


I.  <  *  < 

2  <  ♦  <  2 


(4-11 ) 


If  we  choose  random  numbers  between  0  and  1  and  assign  them  to  s,  the 
corresponding  values  of  $  will  have  the  proper  distribution.  For  SHAD1 , 
N  random  angles  would  be  generated  in  this  manner  and  paired  with  the 
randomly  generated  sighting  times  for  each  orbital  revolution. 

4.3.3  Real  Star  Field  ( SHAD 2 ) 


4. 3.3.1  Star  Catalog 

The  Yale  University  Catalog  of  Bright  Stars  lists  9091  stars  according 
to  position,  type,  and  visual  magnitude.  From  this  catalog,  a  new  catalog 
of  362  stars  was  created  which  contained  those  stars  having  a  brightness  of 
at  least  fourth  magnitude  as  seen  by  an  S-20  star  sensor. 

4. 3. 3. 2  Estimation  of  Sighting  Times 

To  estimate  the  time  when  a  satellite  can  expect  to  see  a  star  rise  or 
set,  the  star's  position  was  first  converted  to  the  orbital  coordinate 
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system  (ip,  _i£,  i^)  shown  in  Figure  4-4.  The  catalogued  star's 
position  as  a  unit  vector  jjg  in  ECI  coordinates  can  be  converted  to 
orbital  coordinates  by  the  transformation 


(4-12) 


where  ip,  _ig  and  are  unit  vectors  defining  the  directions  of  the 
orbital  coordinate  axes  in  ECI  coordinates.  From  this  one  can  define  the 
star's  direction  in  terms  of  an  orbital  right  ascension  a  (measured  from 
the  perigee  direction,  ip)  and  orbital  declination  6  (measured  from  the 
orbital  plane)  as  follows: 


_1  *  * 

a  =  tan  u  /u 
s  /  s 
P  e 


(4-13) 
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For  a  particular  orbit,  the  horizon  elevation  angle  0  is  given  by 


sin  0 


(re  +  byK1  +  6  COS  f) 

P 


(4-15) 


where 

re  =  mean  Earth  radius 
hv  *  vacuum  altitude  of  starlight 
p  =  semi-latus  rectum  of  orbit 
e  =  orbital  eccentricity 

f  =  true  anomaly  (satellite  angular  position  with  respect  to  perigee ) 


To  place  the  Earth  directly  between  the  satellite  and  the  star,  the 
satellite  should  be  180  degree  from  the  projection  of  the  star  direction 
onto  the  orbital  plane,  i.e. 

f  =  a  +  tt  (4-16) 

Substitution  of  this  value  of  f  into  Equation  4-15  yields  a  value  for  0 
which  if  greater  than  the  star's  orbital  declination  5  implies  that  the 
star  will  rise  and  set  at  the  vacuum  altitude  hv  (see  Figure  4-5).  The 
angle  6  shown  in  Figure  4-5  describes  the  angular  half-width  of  the  chord 
formed  by  the  passage  of  a  star  behind  the  Earth.  By  spherical 


trigonometry 


GRAZING  STAR 
(SETTING) 


Figure  4-6.  Orbital  positions  for  observing  a  grazing  star 


and  for  a  setting  star 

fs  =  a  +  it  -  A  (4-20) 

The  satellite  angular  position  where  the  grazing  star  is  observed  may  be 
solved  recursively  using  Equations  4-18  and  4-19  or  4-20.  From  this  the 
time  for  observing  the  star  is  calculated. 


4. 3. 3. 3  Effects  of  Orbital  Precession 

The  estimate  of  satellite  sighting  times  developed  above  presumes  a 
stationary  Kepler ian  orbit  for  at  least  one  revolution.  One  can 
approximate  the  precise  sighting  times  for  a  perturbed  trajectory  by 
considering  only  the  effects  of  the  zonal  harmonic  term  J2«  Since  J2 
dominates  all  other  effects  for  low  earth  orbits,  this  corrected  sighting 
time  should  assure  that  when  the  indicated  star  is  sighted  from  a  perturbed 
trajectory  the  vacuum  altitude  obtained  should  be  very  close  to  the  one 
desired. 

The  presence  of  J2  causes  the  orbit  to  precess  about  the  Earth's 
polar  axis.  The  amount  of  this  precession  during  a  single  orbit  is  given 
by 

A  ft  =  -3ir  J2  (r  /p)  cos  i  (4-21) 

where 

J2  =  gravitational  potential  term 
req  =  equatorial  radius  of  the  Earth 
p  =  semi-latus  rectum  of  the  orbit 
i  =  orbital  inclination 

To  correct  for  orbital  precession  after  each  orbit,  the  orbital 
coordinate  unit  vectors  Up,  i^,  ijj)  used  in  Equation  4-12  were 


rotated  about  the  Earth's  polar  axis  by  the  amount  given  in  Equation  4-21 


SECTION  5 


NAVIGATION  PERFORMANCE  RESULTS 


5.1  Introduction 

Performance  results  were  generated  for  three  different  satellite  orbits 
which  are  considered  to  be  fairly  representative  of  those  used  today.  All 
of  the  orbits  were  circular  and  had  the  following  distinguishing 
characteristics: 


Type  Orbit 

Altitude (km) 

Inclination (deg ) 

Low  Orbit 

915 

65 

GPS 

20,242 

63 

Geosynchronous 

35,866 

0 

The  sensitivity  of  performance  to  variation  in  measurement  error  and 
number  of  star  sightings  is  presented  for  all  of  the  above  orbits. 

However,  the  effect  of  variation  in  other  parameters  is  shown  only  for  the 
low  orbit  and  can  be  considered  to  be  somewhat  similar  for  the  other 
orbits. 

In  most  cases,  only  "setting"  stars  were  used  for  navigation  (i.e. 
stars  in  the  aft  direction  of  satellite  motion).  Such  stars  can  be 
acquired  before  they  undergo  refraction  which  is  much  easier  than  acquiring 
stars  coming  out  of  refraction.  However,  since  there  are  some  performance 

f 

benefits  in  using  both  setting  and  rising  stars,  data  was  also  generated 
for  this  situation  so  that  one  can  weigh  those  benefits  against  the 
additional  difficulty  in  system  design. 
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All  star  sightings  for  position  update  occurred  when  the  ray  grazing 
height  was  25  kilometers.  As  previously  indicated  in  Section  2.4.5,  each 
star /horizon  measurement  was  treated  as  an  equivalent  vacuum  tangent  height 
measurement  so  that  the  errors  in  atmospheric  density  modeling,  and  in 
dispersion  or  refraction  measurement,  could  be  expressed  in  common  units  of 
meters  without  being  concerned  with  the  relative  contributions  of  each. 

Most  of  the  performance  results  were  generated  with  a  baseline  measurement 
error  of  70  meters  (one  sigma)  in  vacuum  tangent  height,  although  some  data 
was  generated  for  other  values. 

The  effect  of  a  bias-type  measurement  error  was  also  included  in  the 
study  and  was  treated  as  an  additional  state  parameter  to  be  updated  by  the 
Kalman  filter.  Again,  no  strong  distinction  was  made  between  the  relative 
bias  error  contributions  made  by  errors  in  atmospheric  modeling  and  in 
dispersion  or  refraction  measurement.  However,  most  individuals  feel  that 
the  major  contributor  will  be  atmospheric  density  modeling,  and  that 
systematic  bias-like  errors  in  atmospheric  modeling  will  occur  over 
temperate  and  polar  regions.  Most  of  the  performance  results  are  for  the 
case  where  a  bias  error  is  present  for  all  observations  outside  +_  30 
degrees  latitude.  At  the  start  of  each  navigation  run,  the  uncertainty  in 
the  bias  estimate  was  generally  assumed  to  be  100  meters  (one  sigma). 
Preliminary  studies  indicated  that  larger  initial  bias  uncertainties  (up  to 
500  meters)  had  little  effect  on  the  final  results. 

It  should  be  noted  that  the  present  treatment  of  a  systematic 
atmospheric  density  error  as  only  a  bias  type  error  is  probably  not  the 
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best  approach  since  one  would  expect  systematic  changes  in  the  error  as  a 
function  of  season  and  geographical  location.  However,  for  the  purposes  of 
this  study,  it  was  felt  that  an  equivalent  bias  in  vacuum  tangdnt  height 
would  be  sufficient. 

Three  different  measurement  bias  options  were  employed  in  the  study  and 
are  as  follows: 

Option  1  -  All  measurements  have  no  bias . 

Option  2  -  All  measurements  have  the  same  bias  which  is  to  be  estimated. 

Option  3  -  All  measurements  within  a  specified  latitude  band . (+  ^ max ) 

have  no  bias  while  those  outside  this  band  have  a  common  bias 
which  is  to  be  estimated. 

» 

5.2  Low  Earth  Orbit 

5.2.1  Introduction 

The  low-earth  orbit  selected  for  study  was  a  circular  orbit  of  915 
kilometers  (497  nautical  miles)  altitude  and  65  degrees  inclination. 
Performance  results  were  generated  using  programs  SHAD1  and  SHAD2  where 
SHAD1  propagated  the  covariances  along  a  Keplerian  orbit  and  SHAD2  included 
the  effect  of  orbital  precession  due  to  the  earth's  zonal  harmonic  term 
J2. 

All  navigation  runs  were  for  a  time  of  flight  of  ten  hours  (*  5.8 

orbits).  At  the  start  of  each  run  the  uncertainty  of  eaci.  component  of 

« 

satellite  position  and  velocity  was  usually  1000  meters  and  1  meter  per 


second,  respectively,  although  other  values  were  sometimes  used  with  very 
little  effect  on  navigation  performance  after  the  first  orbit  of  operation. 

In  addition,  these  uncertainties  were  initially  assumed  to  be  uncorrelated 
which  is  seldom  the  case  in  real  practice  and,  when  included,  usually  reduces 
the  early  transients  in  the  state  estimates. 

The  navigation  runs  performed  with  SHAD1  used  a  time  integration  step  size 
of  30  seconds.  Using  the  program  GTDS  in  the  manner  indicated  in  Section 
4.2.3,  a  value  of  2x1 kilometers/second^  was  obtained  for  the  non-zero 
terms  of  the  process  noise  spectral  density  matrix(Q).  For  a  time  step  of  30 
seconds,  the  corresponding  values  of  the  non-zero  terms  of  the  process  noise 
covariance  matrix  (Q^)  were  6x1 0-1^  kilometers^/second^. 

A  typical  example  of  the  performance  obtained  during  a  navigation  run  is 
show  in  Figure  5-1  where  the  time  histories  are  given  for  the  standard 
deviations  of  the  tangential,  normal  and  radial  components  of  satellite 
position  error.  The  results  are  for  the  case  where  sightings  were  made  on  40 
different  setting  stars  per  orbit  and  the  stars  were  restricted  to  be  within 
45  degrees  of  the  orbital  plane  (i.e.  the  observations  were  restricted  to  an 
azimuth  look -angle  (9)  of  180  45  degrees  about  the  local  vertical  where  zero 

azimuth  is  defined  to  be  in  the  direction  of  satellite  motion).  During  the 
early  part  of  the  first  orbit  it  is  seen  in  Figure  5-1  that  large  position 
error  transients  occur  because  of  the  initial  (and  uncorrelated)  position  and 
velocity  uncertainties  selected  for  the  navigation  run.  However,  as  the 
navigation  filter  continues  to  process  measurements,  the  errors  are  brought 
under  control  during  the  latter  portion  of  the  first  orbit  and  start 
approaching  steady-state  levels.  Note  that  at  least  one  of  the  error 
components  is  oscillatory  in  nature  which  is  actually  the  case  for  all  three 


components  and  is  inherent  in  the  dynamics  of  this  problem.  Since  these 


components  happen  to  oscillate  asynchronously,  it  would  be  ambiguous  to 
quantify  the  navigator's  performance  at  a  specific  time  such  as  at  the  end 
of  the  navigation  run.  A  more  consistent  and  representative  performance 
indication  is  obtained  by  performing  a  linear  least-squares  fit  of  the  data 
over  the  last  orbital  period  and  extrapolating  the  resulting  line  to  the 
final  time.  For  the  case  in  Figure  5-1,  the  average  estimates  of  the 
tangential,  normal  and  radial  components  at  the  final  time  are  39,  43  and 
11  meters,  respectively. 


Number  of  Orbits 


Figure  5-1.  Time  history  of  navigation  position  error  standard  deviations 
(low-earth  orbit) 


5.2.2  SHAD1  Performance  Results 


5 . 2 . 2 . 1  Bias  Option  1  (No  Bias)- 

Before  presenting  data  for  those  cases  containing  measurement  bias  it 
is  felt  that  some  indication  should  be  provided  of  the  performance  when 
there  is  no  bias,  and  no  attempt  is  made  to  estimate  a  bias  along  with 
satellite  position  and  velocity.  For  this  option  the  effect  of  increasing 
the  number  of  star  sightings  per  orbit  was  analyzed  for  two  different 
sighting-azimuth  restrictions. 

Figure  5-2  shows  the  effect  of  increasing  the  number  of  star  sightings 
per  orbit  when  only  stars  within  the  azimuth  sector  of  180  jf  45  degrees  are 
used.  In  addition,  this  figure  provides  some  indication  of  the  effect  of 
restricting  the  observations  to  certain  latitudinal  regions  flf  the  world. 
The  data  on  the  left  of  the  vertical  dashed  line  was  obtained  for 
observations  only  within  _+  30°  latitude,  while  that  on  the  right  was  for 
observations  within  _+  67.5°  latitude.  Note  that  most  of  the  performance 
improvement  is.  achieved  after  making  30  to  40  star  sightings  per  orbit. 

Also  note  that  for  30  or  more  star  sightings  per  orbit  there  is  little  to 
be  gained  by  going  from  a  latitude  restriction  of  _+  30°  to  +  67.5°.  This 
result  tends  to  support  the  idea  of  using  only  observations  within  some 
latitude  band  about  the  equator  where  there  is  less  random  and  systematic 
variation  in  stratospheric  density. 

Figure  5-3  shows  the  sensitivity  of  performance  to  number  of  star 
sightings  per  orbit  when  sightings  are  made  in  both  the  180  +,  45°  and  +_  45° 
azimuth  sectors  (i.e.  both  setting  and  rising  stars  are  used).  It  is  seen, 
by  comparing  Figure  5-3  with  Figure  5-2,  that  little  improvement  is 
achieved  by  including  the  _+  45°  azimuth  sector  except  for  the  normal  compo¬ 
nent  of  position  error  when  the  number  of  star  sightings  per  orbit  is 


small 
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POSITION  UNCERTAINTY  (METERS) 


A  fixed  bias  error  in  all  of  the  observations  can  be  due  to  one  or  more 
'of  the  following: 

1)  Fixed  dispersion/refraction  measurement  error. 

2)  Fixed  time  reference  error. 

3)  Incorrect  knowledge  of  Earth's  radius. 

4)  Fixed  atmospheric  density  error  for  all  observations. 

The  first  two  error  sources  can  be  minimized  by  proper  system  design  and 
the  third  is  considered  to  be  minor.  However,  the  fourth  error  source  is 

m 

considered  to  be  the  most  important  and  likely  source  because  of  present 
limitations  in  man's  knowledge  of  atmospheric  behavior  with  season, 
geographic  location,  etc.  Meteorological  observations  of  the  stratosphere 
near  the  equator  indicate  a  fairly  well-behaved  and  predictable  density 
with  little  chance  of  a  meaningful  bias  error.  As  one  departs  from  the 
equator  by  more  than  20  or  30  degrees  latitude,  there  are  larger  systematic 
changes  which  have  some  meaningful  uncertainty.  For  this  reason  Option  3, 
presented  later,  is  felt  to  be  more  representative  of  the  situation  than 
Option  2.  However,  as  a  matter  of  interest,  some  performance  data  was 
generated  for  Option  2. 

Figure  5-4  shows  the  results  with  Option  2  when  star  sightings  are  made 
within  the  azimuth  sector  of  1 80  +_  45  degrees.  There  was  no  latitude 
restriction  on  the  star  sightings.  Note  that  the  performance  in  estimating 
the  normal  and  radial  components  is  essentially  the  same  as  Option  1  (see 
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Figure  5-2).  However,  the  tangential  position  error  of  80  to  110  meters  is 
twice  as  large  as  that  for  Option  1 .  The  reason  why  the  tangential 
performance  is  not  as  good  is  due  to  the  difficulty  of  the  navigation 
filter  in  distinguishing  between  a  tangential  position  error  and  a 
measurement  bias  error  when  using  only  star  sightings  in  the  aft 
direction.  A  larger  azimuth  sector  in  the  aft  direction  and/or  more  orbits 
of  data  processing  would  undoubtedly  result  in  further  improvement  of  the 
tangential  component,  however,  a  more  dramatic  improvement  may  be  achieved 
by  also  using  forward  star  sightings.  Use  of  both  forward  and  aft  star 
sightings  enables  the  navigation  filter  to  easily  distinguish  between  the 
tangential  position  error  and  a  measurement  bias  error.  This  is  clearly 
demonstrated  in  Figure  5-5  where  star  sightings  were  made  in  the  180  +_  45° 
and  _+  45°  azimuth  sectors.  The  overall  prerformance  is  equivalent  to  that 
obtained  with  no  measurement  bias  and  demonstrates  the  ability  of  the 
filter  to  handle  a  bias. 

5. 2. 2. 3  Bias  Option  3 

Option  3  assumes  that  all  measurements  within  a  specified  latitude  band 
(+^max)  have  no  bias  while  those  outside  have  a  common  bias  which  is  to 
be  estimated.  As  previously  indicated,  this  option  is  considered  to  be  the 
most  representative  of  the  three  options  if  the  bias  is  associated  with 
atmospheric  density  modeling. 

Figure  5-6  shows  the  performance  with  Option  3  when  the  azimuth  sector 
is  180  _+  45°  and  a  bias  exist  only  outside  _+  30°  latitude.  Note  that  the 
performance  is  similar  to  that  obtained  with  no  bias  anywhere  (Figure  5-2) 
even  though  the  star  sightings  are  restricted  to  the  180  _+  45°  azimuth 
sector.  It  would  appear  that  no  sightings  are  required  on  rising  stars 
(i.e.  stars  in  the  forward  azimuth  sector)  if  there  is  no  atmospheric  bias 
error  within  +  30°  latitude. 
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Figure  5-5.  Sensitivity  of  performance  to  number  of  star  sightings  when 
all  measurements  are  biased  but  are  made  fore  and  aft 
(low-earth  orbit). 
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Shown  in  Figure  5-7  is  the  effect  of  changing  the  azimuth  sector  size 
for  star  sightings  in  the  aft  direction.  The  results  are  for  40  star 
sightings  per  orbit  and  no  bias  within  _+  30°  latitude.  Note  that  only  the 
normal  component  is  significantly  improved  by  using  a  sector  size  larger 
than  180  _+  45  degrees.  This  is  as  expected  since  star  sightings  at  larger 
azimuth  angles  with  respect  to  the  orbital  plane  provide  more  normal 
positional  information.  However,  it  should  also  be  noted  that  stars  at 
large  azimuth  angles  require  more  time  to  enter  the  atmosphere  than  those 
near  the  orbital  plane. 

A  baseline  measurement  error  of  70  meters  (one  sigma)  was  used  to 
generate  all  of  the  prior  performance  data.  Figure  5-8  shows  the  effect  of 
varying  this  error  for  the  case  of  40  star  sightings  per  orbit,  no  bias 
within  _+  30°  latitude,  and  an  azimuth  sector  of  180  _+  45°.  It  is  seen  that 
the  effects  are  fairly  linear,  with  the  radial  error  remaining  small  and 
the  other  two  components  being  about  half  as  large  as  the  measurement 
error. 

5.2.3  SHAD2  Performance  Results 

The  program  SHAD2  provides  the  capability  to  model  satellite  orbits 
more  accurately  and  to  base  star  observations  upon  real  star  positions. 
Since  the  zonal  harmonic  coefficient  J2  dominates  all  other  disturbing 
acceleration  coefficients  in  low-earth  orbit,  it  was  decided  to  study  the 
effect  of  this  term  alone  on  navigation  preformance.  The  star  positions 
were  taken  from  the  Yale  Bright  Star  Catalog.  , 

Figure  5-9  shows  the  effect  of  varying  the  number  of  star  sightings  per 
orbit  with  SHAD 2  for  the  case  where  the  other  conditions  were  the  same  as 
those  used  to  generate  the  SHAD1  results  of  Figure  5-6  (i.e.  star  sightings 


Figure  5-7.  Sensitivity  of  performance  to  azimuth  sector  size 
(low-earth  orbit) . 
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Figure  5-9.  Sensitivity  of  performance  to  number  of  star  sightinqs  when 
using  SHAD2  program  (low-earth  orbit) . 


restricted  to  an  azimuth  sector  of  180  _+  45°  and  a  measurement  bias  only 
outside  +_  30°  latitude).  It  is  seen  that  the  results  are  about  the  same  as 
those  in  Figure  5-6,  thus  providing  some  validity  to  the  results  of  SHAD1 . 
These  results  were  also  confirmed  by  M.  Slutsky  using  the  Kalman  filtering 
options  of  GTDS .  As  a  matter  of  interest,  the  left-most  set  of  data  points  in 
Figure  5-9  represent  the  performance  near  the  end  of  the  time  history  in 
Figure  5-1 . 

5.3  GPS  Orbit 

The  second  orbit  analyzed  in  the  study  is  similar  to  that  of  the  Global 
Positioning  System  (GPS)  and  had  an  altitude  of  20,242  km,  orbital  inclination 
of  63°,  and  period  of  12  hours. 

All  of  the  performance  data  generated  for  this  orbit  was  obtained  using 
program  SHAD1 .  The  covariance  matrix  was  propagated  at  150  second  intervals 
for  100  hours  (8.3  orbits).  Using  the  program  GTDS  in  the  manner  indicated  in 
Section  4.2.3,  a  value  of  1 .4  x  10-17  kilometers2/second2  was  obtained 
for  the  non -zero  terms  of  the  process  noise  covariance  matrix  ( ) . 

A  typical  example  of  the  performance  during  a  navigation  run  is  shown  in 
the  time  history  of  Figure  5-10.  The  results  are  for  the  case  of  26  star 
sightings  per  orbit  with  no  sightings  outside  an  azimuth  sector  of  180  ±  90°, 
no  observations  outside  +_  30°  latitude  and  no  bias  error.  It  is  seen  that  the 
initial  transients  settle  out  after  one  orbit  of  data  processing. 

Figure  5-11  shows  the  effect  of  changing  the  number  of  star  sightings  per 
orbit  when  all  other  conditions  remain  the  same  as  those  of  Figure  5-10.  The 
right-most  data  set  of  Figure  5-11  represents  the  final  average  values  of  the 
time  histories  in  Figure  5-10. 
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Figure  5-10.  Time  history  of  navigator  position  error  standard  deviations 
(GPS  orbit) . 


Figure  5-12  shows  the  effect  of  changing  the  number  of  star  sightings 
per  orbit  when  all  other  conditions  are  the  same  as  those  of  Figure  5-11 
except  that  observations  are  also  permitted  outside  _+  30°  latitude  but 
contain  a  bias  error  which  is  to  be  estimated.  No  latitude  restrictions 
permit  observations  to  be  made  throughout  the  orbit,  however,  it  is  seen 
that,  for  20  or  more  star  sightings  per  orbit,  the  performance  is  only 
slightly  better  than  that  of  Figure  5-11  where  all  of  the  observations  were 
inside  ±  30°  latitude. 
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Figure  5-11.  Sensitivity  of  performance  to  number  of  sightings  when  all 
measurements  are  within  +  30°  latitude  (GPS  orbit) . 
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The  sensitivity  of  performance  to  measurement  error  is  shown  in  Figure 


5-13  for  the  case  of  10  star  sightings  per  orbit  and  for  an  azimuth  sector 
of  180°  ±  45®.  Note  that  the  normal  component  of  position  uncertainty  is 
now  larger  than  the  tangential  component  primarily  because  the  azimuth 
sector  is  smaller  than  that  used  for  the  results  in  Figures  5-11  and  5-12. 

5.4  Geosynchronous  Orbit 

The  third  orbit  analyzed  in  the  study  was  a  geosynchronous  orbit  which 
lies  within  the  equatorial  plane  and  has  an  altitude  of  35,866  km  and 
orbital  period  of  24  hours. 

All  of  the  performance  data  generated  for  this  orbit  was  obtained  with 
the  SHAD1  program.  The  covariance  matrix  was  propagated  at  300  second 

m 

intervals  for  200  hours.  Using  the  program  GTDS  in  the  manner  indicated  in 
Section  4.2.3,  a  value  of  1 .5  x  10-20  kilometers2/second2  was  obtain¬ 
ed  for  the  non-zero  terms  of  the  process  noise  covariance  matrix  (Q^). 

All  of  the  performance  results  presented  for  the  geosynchronous  orbit 
are  for  Bias  Option  3  (see  Section  5. 2. 2. 3)  where  all  observations  within  a 
specified  latitude  band  (±  30°)  h;?ve  no  bias  while  those  outside  have  a 
common  bias  which  is  to  be  estimated. 

Figure  5-14  shows  the  navigator  performance  as  a  function  of  the  number 
of  star  sightings  per  orbit.  Since  the  satellite  velocity  is  much  slower 
than  at  low-earth  orbit,  the  amount  of  divergence  in  the  error  covariance 
matrix  between  measurements  is  also  much  slower.  Thus,  the  reduced  number 
of  grazing  stars  seen  from  this  altitude  does  not  sacrifice  performance 
compared  to  low-earth  orbit. 
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The  sensitivity  of  performance  to  measurement  error  is  shown  in  Figure 


5-15  for  the  case  of  10  star  sightings  per  orbit.  Here  again,  it  is  seen 
that  the  sensitivities  are  similar  to  those  obtained  for  the  GPS  orbit, 
with  the  normal  component  of  position  uncertainty  being  about  twice  as 
large  as  the  tangential  component. 
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Figure  5-15.  Sensitivity  of  performance  to  measurement  error 
(geosynchronous  orbit) . 


APPENDIX  A 


Derivation  of  ha(R,  pQ)  for  an  Exponential  Atmosphere 

The  density  of  am  exponential  atmosphere  is  completely  defined  by  the 
following: 


( A— 1  ) 


where  p  is  the  density  at  altitude  h,  pQ  is  the  density  specified  at  some 
arbitrary  altitude  hQ/  and  H  is  the  specified  constant  for  density  scale 
height.  According  to  Equation  2-6  the  approximate  refraction  for  this 
atmosphere  is  as  follows: 

»-»<»■>,[*  fry-h,1]T 


where  hg  is  the  grazing  height  of  the  refracted  ray  and  pg  is  the 
density  at  that  height.  Since  re  >>  hg  the  refraction  may  be  expressed 
as  follows: 

B  2  (A'3’ 


Using  Equation  A-1  the  density  pg  can  be  expressed  as 


which  when  substituted  for  pg  in  Equation  A-3  yields 


R  *  MX  )p 


( A-5 ) 


Inverting  the  above  expression  gives 
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According  to  Equation  2-18  the  relationship  between  the  apparent  height 
(ha)  and  grazing  height  (hg)  for  a  spherically  stratified  atmosphere  is 


['  »,] 


h  »  1  +  k  (X  )  p  h  +  k (x  )p  r 

a  g  g  g  e 


(A-7 ) 


Since  the  term  MX  )pg  is  never  any  larger  than  about  2  x  10“^  for 
altitudes  above  20  kilometers,  it  may  be  neglected  in  the  first  expression 
on  the  right  of  Equation  A-7  so  that 


h  •>  h  +MX)p  r 

a  g  g  e 


(A-8) 


Substituting  for  pg  from  Equation  A-3  and  for  hg  from  Equation  A-6 
gives  . 


h  (R,  p  )  -  h  -  H  ln(R)  +  H  In 
a  o 


MX  )p 


fi-!s 

>L  H 


+  R 


H  r 


L  2u 


(A-9 ) 


-  116 


a 


9 

a 


>3 


$ 


lyf 

'a 


8 


LIST  OF  REFERENCES 


Paulson,  D.C.,  "Autonomous  Satellite  Navigation  From  Strapdown 
Landmark  Measurements,"  Proc.  3rd  Symp.  on  Nonlinear  Estimation  Theory 
(San  Diego,  CA,  September  1972),  p.  167-182. 


Lowrie,  J.W.,  "Autonomous  Navigation  Systems  Technology  Assessment,” 
AIAA  Report  No.  79-0056,  January  1979. 


Fang,  B.T.,  "Satellite-To-Satellite  Tracking  Orbit  Determination," 
Journal  of  Guidance  and  Control,  Vol.  2,  No.  1,  Jan. -Feb,  1979, 
p.  57-64. 


Marshall,  M.H.,  "Goals  for  Air  Force  Autonomous  Spacecraft,"  Jet 
Propulsion  Laboratory,  Report  No.  SD-TR-81-72,  1981. 


Pease,  G.E.,  and  H.T.  Hendrickson,  "Autonomous  Navigation  Accuracy 
Using  Simulated  Horizon  Sensor  and  Sun  Sensor  Observations,"  Aerospace 
Corporation,  Los  Angeles,  CA,  (Available  NTIS  SAP:HC  A20/MF  A01), 
October  1980. 


Hoag,  D.G.,  "The  History  of  Apollo  On-Board  Guidance,  Navigation,  and 
Control,"  (presented  at  the  International  Space  Hall  of  Fame 
Dedication  Conference,  October  1976),  The  Charles  Stark  Draper 
Laboratory,  Inc.,  Cambridge,  MA,  Report  No.  CSDL-P-357,  1976. 


Lillestrand,  R.L.,  and  J.E.  Carroll,  "Horizon-Based  Satellite 
Navigation  Systems,"  IEEE  Transaction  on  Aerospace  and  Navigational 
Electronics,  September  1963,  p.  247-270. 


Barnes,  F.,  "Earth  Horizon  Sensing  by  Stellar  Refraction  Measurement," 
The  Charles  Stark  Draper  Laboratory,  Inc.,  Cambridge,  MA,  Intralab 
Memorandum  No.  AGS-25-75,  16  October  1975. 


Elden,  B.,  "The  Dispersion  of  Standard  Air,"  Journal  of  pptical 
Society  of  America,  Vo.  43,  No.  5,  May  1953,  p.  339-344. 


U.S.  Standard  Atmosphere,  1976,  NOAA  -  S/T  76-1562,  NOAA,  NASA,  U.S 
Air  Force,  Washington,  D.C.,  October  1976. 


11.  Barnes,  F.,  "Starlight  Refraction  by  Ray-Tracing  thru  Modeled  Earth 
Atmospheres  Using  the  "RFRAC3,"  Program,  "The  Charles  Stark.  Draper 
Laboratory,  Inc.,  Cambridge,  MA.  Intralab  Memorandum  No.  AGS-19-81, 
2  November  1 981 . 


12.  Cole,  A.E.,  and  A.J.  Kantor,  "Air  Force  Reference  Atmospheres,"  Air 
Force  Geophysics  Laboratory,  Hanscom  AFB,  MA,  Report  No. 
AFGL-TR-78-0051 ,  28  February  1978. 


13.  Cole,  A.E.,  A.J.  Kantor,  and  C.R.  Philbrick,  "Kwajalein  Reference 

Atmospheres,  1978,"  Air  Force  Geophysics  Laboratory,  Hanscom  AFB,  MA, 
Report  No.  AFGL-TR-79-0241 ,  24  September  1978. 


14.  Gelb,  A.,  (ed.),  "Applied  Optimal  Estimation,"  MIT  Press;  Cambridge, 
MA,  1974. 


15.  Goodyear,  W.H.,  "Completely  General  Closed-Form  Solution  for 
Coordinates  and  Partial  Derivatives  of  the  Two-Body  Problem," 
Astronomical  Journal,  Vol.  70,  No.  3,  April  1965,  pp.  189-192, 


16.  Gai,  E.,  and  R.C.  Hutchinson,  "On-Board  Navigation  Equations,”  The 
Charles  Stark  Draper  Laboratory,  Inc.,  Cambridge,  MA,  Report  No. 
CSDL-C-5360,  November  1980. 


17.  Schweppe,  F.C.,  "Uncertain  Dynamic  Systems,"  Prentice-Hall,  Englewood 
Cliffs,  New  Jersey,  1973. 


18.  Computer  Sciences  Corporation  and  Systems  Development  and  Analysis 
Branch  (GSFC),  joint  ed.,  "Goddard  Trajectory  Determination  System: 
User's  Guide,  "  July  1978. 


19.  Battin,  R.H.,  "Astronomical  Guidance, "  McGraw-Hill,  New  York,  NY, 

1964. 


